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ABSTRACT 


This report presents new mass transport elements for the next generation of the NIST IAQ Model 
that may be used to model a) homogenous (bulk-air) chemistry within well-mixed chamber, b) 
aerosol mass transport within well-mixed chambers and fractional particle filtration in building 
filtration devices,and c) heterogeneous (surface-related) physical processes and chemical 
transformations including those governing the behavior of gas-phase air cleaning devices. In an 
effort to maintain rigor, generality, and flexibility, each transport process is formulated in terms of 
the elemental mass transport steps that together govern the overall process. In this way, the more 
complex processes may be represented aS component equations that are assembled from 
fundamental element equations. 


The element/component assembly method, upon which the NIST IAQ Model is based, provides a 
general and modular approach to the formulation of systems of equation governing the mass and 
air transport in buildings to effect indoor air quality analysis. In this approach, the solution of the 
system equations is a computationally distinct task that may be achieved using a variety of 
numerical methods. The third chapter of this report discusses numerical and computational 
strategies for the solution of the system equations that are compatible with both the existing and 
proposed new mass transport elements and presents candidate strategies that appear to be most 
promising. 


Finally the fourth chapter of this report considers user interface strategies to implement the 
proposed new mass transport elements and components. 


Keywords: contaminant dispersal, filtration, indoor air quality, mass transport, modeling, 
ventilation 
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1. INTRODUCTION 


The NIST IAQ Model has now undergone four generations of development marked by the release 
of the programs CONTAM86, CONTAM87, CONTAM88, and, most recently, CONTAM93/94. 
CONTAM86 supported single species multizone gas-phase contaminant dispersal analysis, 
CONTAM87 extended the NIST IAQ Model to account for multiple reactive gas-phase 
contaminants and the possibility of convection-diffusion mass transport in 1D flow regimes (e.g., 
ducts), CONTAM88 integrated multizone airflow analysis with the contaminant dispersal analysis 
of the earlier generations and CONTAM93/94 was developed to provide users with a graphic users 
interface, more complete libraries of mass transport and airflow elements and better numerics. 


During the past eight years of model development, the IAQ research community has not stood still. 
More complete mathematical models have emerged to account for a) homogeneous (bulk-air) and 
heterogeneous (surface-related) physical and chemical transformations, b) aerosol/particulate 
transport, c) pollutant source emissions, and d) sorption and particle filtration devices. 
Consequently, the time is ripe to begin the development of new mass transport elements and 
components for the next generation NIST IAQ Model based on these more complete models. 


This report is directed toward the development of these needed mass transport elements and is 
organized into three chapters — formulation of element equations, identification of strategies for 
solution of system equations and development of user interface strategies for computational 
implementation. 


In preparation for this report, the literature was reviewed to assemble a) available mathematical 
“models for the transport processes enumerated above and b) user interface strategies used to 
implement these models computationally. From these reviews the most promising models for 
homogeneous chemistry, aerosol transport, and filtration were selected. For source emission 
modeling, on the other hand, a new approach was taken. Elemental transport relations were 
collected for the fundamental transport steps that may be expected to govern the actual physical 
mechanisms governing emission for a variety of possible sources. These elemental relations were 
then used to assemble a series of examples of physically-based emission models without any 


attempt to be comprehensive. Empirically-based emission models were not considered. 


Few of the available models are presented in the necessary element equation form most suitable for 
computational implementation within the CONTAM element assembly framework. Therefore, the 
available models were (re)formulated as element equations. In an effort to maintain rigor, 
| generality, and flexibility, each transport process was reformulated in terms of the elemental mass 
transport steps that together govern the process. In this way, the more complex processes may be 
represented as component equations that are assembled from fundamental element equations. 
Thus, for example, sorption filtration components are developed as assemblages of elements 
modeling boundary layer diffusion, pore diffusion, and sorption. 
5 


The element assembly method, upon which the CONTAM family of program is based, provides a 
general and modular approach to the formulation of systems of equations governing the mass and 
air transport in buildings to effect indoor air quality analysis. In this approach, the solution of the 
system equations is a computationally distinct task that may be achieved using a variety of 
numerical methods. The choice of numerical method used is, however, limited by the 
mathematical character of element equations being assembled. The third chapter of this report 
discusses numerical and computational strategies for the solution of the system equations that are 
compatible with both the existing and proposed new mass transport elements and presents 


candidate strategies that appear to be most promising. 


The methods used to communicate input to a program (1.e., building idealization, representation, 
and program interface conventions) and computed results to the user demand careful consideration. 
The last chapter attempts to develop user interface strategies, based on the review of conventions 
used in related programs, for the implementation of the proposed new mass transport elements. 


Throughout the report, text marked, as this paragraph is marked, by a V in the left margin and a 
dotted line in the right margin provide specific recommendations for the next generation of the 
NIST IAQ Model — CONTAM«xx. 


2. ELEMENT/COMPONENT EQUATIONS 


Fundamental theory and equations derived to be integrable with current multizone indoor air quality 
analysis models — most importantly, the NIST IAQ Model — will be presented in this chapter to add 
mass transport modeling capabilities for: 


e Homogeneous Chemistry 
¢ Aerosol Transport & Fractional Particle Filtration 


e Heterogeneous Processes 
— Elemental Mass Transport Processes 
— Room Sorption Transport 
— Sorption Filtration 
— Fundamentally-Based Air Pollutant Source Models 


User interface requirements and numerical methods to form and/or solve systems of equations 
assembled from the element or component equations presented will be discussed within the 
individual sections of this chapter. Chapters 3 and 4 of this report will then attempt to review these 
individual interface requirements and numerical methods to formulate strategies to efficiently a) 
form and solve systems of equations assembled from these new and existing element or component 
models and b) communicate input data and output results with the users of programs based on 
these new and existing models. 

2.1. Homogeneous Chemistry 


Indoor air quality may be affected by chemical reactions that occur homogeneously within the bulk 
air phase! contained within building rooms or zones. Three broad classes of homogeneous 
chemistry will be distinguished here — gas-phase, aqueous-phase, and aerosol-phase chemistry — 
where aqueous-phase homogeneous chemistry is a special, but important, subclass of aerosol- 
phase chemistry involving chemistry occurring within or on individual particles of a 
homogeneously distributed aerosol. It is presumed, here, that the dynamics of these three classes 
of homogeneous chemistry (i.e., as controlled by both the kinetics of the chemistry in question and 
the mass transport processes that may affect this chemistry) may be modeled using the vast body of 
theory, knowledge, and information developed in the closely related field of (outdoor) 
tropospheric, atmospheric chemistry (1-4). 


This report will limit consideration to the development of new gas-phase homogeneous chemistry 
mass transport elements for the next generation NIST IAQ model leaving aqueous-phase and 
aerosol-phase chemistry to future efforts. Recent research has indicated that gas-phase chemistry 


! bulk air-phase — used in chemical engineering literature to refer to the main or greater part of a volume of air as 
distinguished from that part of the air volume associated with surface layers at the boundaries of a volume of air. 
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may play a significant role in indoor air at times (5-10) and aqueous-phase chemistry is presently 
coming under closer scrutiny (11-13). For this reason, issues relating to aqueous-phase and 


aerosol-phase chemistry will be reviewed to aid the planning of future development efforts. 
2.1.1. Gas-Phase Chemistry 


Fortuitously, the detailed mechanism and associated kinetics of gas-phase atmospheric chemistry 
can, most often, be described in terms of either unimolecular, bimolecular, or termolecular 


elementary reactions and their kinetic rate expressions as: 


ky 
Unimolecular A — products (la) 
GlA\t eee 
a = -ky [A] (1b) 
k 
Bimolecular n,A + npB + ncC + npD + ... (2a) 
1 d[A] _ 1 a[B) _ meepdlCh me atalD) & 
SE ee AAS Ne dix Ini tie ade (2b) 
k 
Termolecular n,A + npB + ncC > npD + ngE + ... (3a) 
nina (A) ee ieal Ble fd Cle _ —_] d[ Dom saley 
On he ATS Vier reste EARNEST Np dtu Ne dtxqual (3b) 


where [A], [B], ... are concentrations conventionally expressed in terms of either molecules- 
species-cm-air™> or ppm (i.e., molecules-species-10®-molecules-air-! or mole-species: 106-mole-air- 
1); n4, np, ... are the number of molecules or moles stoichiometrically involved in each reaction; 
and kj, k2, k3 are first, second, and third order rate constants respectively. 


Given one mole or 6.023 x 1023 molecules occupies a volume (Vs) of 22,400 cm3 at standard 
conditions of temperature (7; = 273 °K) and pressure (Ps = 1 atm) or a volume V at any other state 
of temperature 7 (°K) and pressure p (atm): 


y = ¥(BIE) : 


the conventional units of the rate constants used above — k;, kz, and k3 — and the conversion factor 
between these units may be derived as tabulated below (see (1) Appendix 4.A.1 for details): 


Table 1: Conventional kinetic rate constant units and their conversion. 


Conversion Relation 


y122:400/p) (7/273) x10° 
6.023102 


Conventional Units 


molecule-cm-3 


koppm = 4.40x10'7(F) ko 


k3ppm = 3-23x10°(E) k 


cm®-molecule-2-s-! 


2.1.1.1. Unimolecular Reactions & Photochemical Dissociation 


Unimolecular elementary reactions often involve a reactant that is in an excited or activated state, 
A*. In atmospheric chemistry, photochemical conversion resulting from the absorption of a 
photon, Av, provides an important path to activation: 


Photochemical Activation A+hv —- A* (5) 


Of the photolytically activated reactions possible, photolytic dissociation (i.e., the fragmentation or 
separation of the original species A into component constituents) is most important: 


Photochemical Dissociation A + hv > A* > npB + n(C +... (6) 


The photolytic dissociation of trace air pollutants may be represented as a first order process, 
Equation 1A, with the rate constant, k;, dependent on the integral of the product of the actinic 
irradiance density , I(A) (photons-cm->-sec-!) — a measure of the spectral distribution of radiant 
flux or radiant intensity (photons-cm-2-sec-!) per unit wavelength (cm-!), the absorption cross- 
section of the molecule A, 6,(A, 7) (cm2), and the quantum yield or probability of dissociation 
of A upon absorption of a photon, ,(A, 7) , over the wavelength interval of interest, (Ay u ho) 
(cm), as: 


Photochemical Dissociation ki = f° O,(A, T)o,4(A, T)I(A) dr (7) 
1 


where T is included to account for the temperature dependency of these quantities. Seinfeld notes 
that in the troposphere (Ay s ho) = (280 730 nm) and for practical application Equation 7 is 


approximated by a discrete approximation using wavelength intervals, AX , of from 10 to 20 nm 
as ( (1) page 113): 
Photochemical Dissociation ky = LOA; TOA» TTA) Ar (8) 


where the overbars are used to indicate mean values taken from measured data. Table 4.1 in 
Seinfeld’s text, lists the major photochemical dissociation reactions occurring in the lower 
atmosphere. 

The kinetics of a unimolecular reaction may be represented as a pseudo, “zero-order” process with 


rate kp,, when the reactant is abundant and the rate of reaction is so slow that the reactant’s 
concentration does not change significantly during the course of the reaction: 


Pseudo Zero Order 1 4dlAl Koes + Koeo = k,[A] for [A] =[A],, = constant (9) 


This condition may also be identified as a pseudo-steady state condition and thus the ss subscript 


notation above for the reactant [A],, . 
2.1.1.2. Bimolecular and Termolecular Reactions 


In atmospheric chemistry, the bimolecular and termolecular reactions described above most often 
involve equal numbers of reactant groups — that is, unit stoichiometry. This is due to the fact that 
in the gas-phase these elementary reactions are invariably associated with collisions of individual 
atoms, molecules, radicals?, or ions. Given the probability of simultaneous collision of three 
atoms, molecules, radicals, or ions is likely to be small, termolecular reactions are uncommon and 
are likely to proceed at low rates. That is to say, atmospheric chemistry is dominated by 
(unimolecular and) bimolecular elementary reactions. 


Tables 4.10 and 4.11 in Seinfeld’s text, lists the important reactions of the oxides of nitrogen 
occurring in the lower atmosphere showing that they are dominated by bimolecular kinetics and 
unit stoichiometry of reactants. A closer examination of these tables reveal two generic cases of 
bimolecular kinetics — that presented above in Equations 2 and the special case: 


A +A — products ; 4ctA) = -k,[A]* (10) 


(Note the stoichiometric coefficient of 1/2 here.) By generalization, three generic cases of 
termolecular kinetics may be identified — that presented above in Equations 3 and the following two 


special cases: 
A+A+A —~products ; icl4) = -k3[A]° (11) 
(A + A) + ngB - products ; peel cs ~k3[A]*[B] (12) 


2 Radicals may be thought to be uncharged (i.e., nonionic) molecular fragments that contain an unpaired electron 
and, as a result, tend to be highly reactive. ) 
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Finally, a special class of termolecular reaction should be noted — the generic combination reaction 
involving a third body M of the form: 


Third Body Combination A+B+M —>5AB+M (13) 


The mechanisms responsible for this type of reaction are thought to involve a series of two 
bimolecular elementary reactions (see (1) Appendix 4.A.2 for details) and the rate expression for 
this reaction is often written in the form of a bimolecular rate expression (i.e., Equation 2b) with 
the rate constant dependent on the concentration of the third body at the low-pressure limit , 
independent of the third body concentration at the high-pressure limit , and a blended combination 
between as: 


k»5,[M] ; as[M] > 0 


Third Body Combination poe eon F (14) 
20 + R200 


i) 
i 


k>.. 5 as[M] > 


where F is the so-called broadening factor used to fine-tune the intermediate relation. That is to 
say, at the low pressure limit this generic combination reaction is governed by termolecular kinetics 
while at the high pressure limit the behavior is bimolecular. 
In a similar way the kinetics of a bimolecular reaction may be represented as a pseudo, first-order 
process with rate constant k,,., when one of the reactants remains practically constant (e.g., if the 
reactant is so abundant that its concentration does not change significantly during the course of the 
reaction) — another example of a pseudo steady state: 
Pseudo First Order 
a[A d([B 
Boks = a = —k,..[A] ; kj.=k[B],, for [B] =[B],, = constant (15) 


2.1.1.3. Temperature Dependencies 


For ideal elementary reactions, rate constants are dependent on the absolute temperature at which 
the reaction occurs as defined by the Arrhenius Equation: 


—E 
k= Aex a2] de (16) 


where A is the Arrhenius constant with the same units as k, Eg; is the so-called “activation energy” 
associated with the given reaction, R is the universal gas constant (8.3145 J-mol-!.°K-!), and T is 
the temperature of the reaction (°K). See (14) pages 14 to 16 for a general introduction to the 
Arrhenius relation. 


Several specific Arrhenius relations are included in the tables of Seinfeld’s text mentioned above. 
A close examination of these tables will reveal, however, other temperature dependent — and more 


1] 


complex — relations (e.g., for the important reaction of oxygen with the oxygen radical to produce 
ozone: 0 + O2 + M 3 03 +M). 


2.1.1.4. Other Rate Relations 


Although it is believed that practically all chemical reactions can be represented in terms of 
unimolecular, bimolecular, and/or termolecular elementary mechanisms these mechanisms may not 
be known or researchers may have had to limit consideration to the overall reaction that results 
from the combined elementary steps. As a result, for a few cases only overall or completely 
empirical rate expressions may be available. As might be expected, such rate expressions assume a 


large variety of forms and, therefore, must be considered on a case-by-case basis. 
2.1.1.5. Formulation for IAQ Analysis — Well-Mixed Zone Idealization 


Thus far, all rate expressions have been presented in the conventional format utilized in the field of 
kinetics. For indoor air quality analysis, the conventions are different — we seek to consider the 
dynamic change of zone (mean) air pollutant concentrations C4, Cp, Cc, ..., expressed in terms of 
mass fraction (mass-A-mass-air™!, etc.), within a collection of well-mixed zones. The mass 
fraction of a species B, for example, is related to the concentration expressed in terms of 
molecule-cm-3 as: 


a 2 ee 
Ca od Se yyy) eet l ( mol-B (Same =} 1 ( cm>—air ] 


grams-—air cm3—air }N,\molecules—B mol-B_ }Pair\ grams—air 
or 
Mz 
Cp =AB 17 
Z INapar ( 


where Ng is Avogadro’s constant (6.022 x 1023 molecules-mol-!), Mg is the molecular weight of 
species B (grams-B-mol-B-!), and pair is the density of air within the well-mixed zone (grams- 
air-cm->). For non-trace analysis one must be careful to use the (current value of the) actual 
density of the pure-air/air-pollutant mixture for Pgir. 


For a well-mixed zone containing a volume of air of mass Mgjr the mass rate of production (+) or 
removal (-) of a given air pollutant, say for species B again, Rg (grams-B-sec-!) is related to the 
time rate of change of species concentration [B] as: 


R a M girM pz d[B) 


N Pair dt ee 


Using Equations 17 and 18 one may, then, convert each of the rate expressions presented above in 


conventional kinetics form to well-mixed zone mass generation or removal rate expressions as 
tabulated below: 


12 


Table 2. Well-Mixed Zone Chemical Generation/Removal Rate Expressions 


Zero Order 


First Order 


Second Order 
Reactant A 


Second Order 
Reactant B 


aA | wy sir 
ate. ko 


aAl = —n,k[AIBI 


All = —npk[A][B] 


M iM, 
Ri a k 
é NAG 
K M k 
a NADair 


MGiyM a 
Napa 


N : 
ky APairc 
M4 


Ra = 
Rag = —M4j,k\Ca 


Kig =k, 
R, = sibel se AC, 


M girM a, NaPair 


R, =-n 
- GAN 4 Ogi M, 


NaPair 
—— €,,C 
Moen 


NaPai 
air Mal kyC4 Cz 


Ry = —Mg;,Ko4CaCp 


n M airMp N Pair NaPair 6 


Rz = - 
a ZINGOn, : M, Mp : 


Cp 


N ; 
Reena Arak xC4Cp 


air 
A 
NaPai 
Ky = ng ky 
M4 


Rp = —Mqi-KopCaCp 


Second Order d{C] 


yt Pair NaP 
= n-k.[A][B] = etic, cai el) cape | ee Ge 
Product C dt pee feb NDR 77M Mp 


Nae 
Rc = nM aioM cg gg F2CaCe 
APair ¢ 
AM 


N 
K5c = = Heh esy 


Rc = Mgi-KocCaCp 


Third Order 
d{A] — —n,k3[A][B][C] Ry =-—n, MairM a, NaPair NaPair Na Pair @ 7;CaCe 


Ie 
Reactant A ” NaPair Mz Msg Mc 


R = -Nn M : k3C,CpCc 
A A air MM 
(NaPair) 
K2, = ny———k 
BA A MM c 3 


R, se a ss aCaCpC 


Third Order d(C] 


——— M,;Mc, N ar Aas a Nipe 
Product C dt = nck3[A][B][C] Rc = nek AP AP AP. NaPairc CC, 


Niogee Mae M, Msg Me 
(NaPair) 
Rc = ncM Mom ks G.ee 
al M,Mpz 


(N A Parl 


K3c = ncMc MM j 


k3 


Ro = MgirK3cCaCpCc 


Other cases follow, by example, from these generic possibilities. 
2.1.1.6. Example Application — Single Well-Mixed Zone with Bimolecular Chemistry 


It is useful to place these equations in the full context that they will appear in indoor air quality 
analysis. To do so, consider a single well-mixed zone containing two reactants A and B, their 
products C and D, and a volume of air of mass Mgjr. For trace dispersal analysis (and no 
significant generation of air within the zone) mass conservation demands that the mass flow rate of 
air exiting the zone be equal to that flowing in, wgjr. Given these conditions, mass balances for 
each of the four pollutant involved may be directly formed as: 
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C4 Ra C4 Cao Gy E, 


Cp Rp ile a \uinanet JAC Bo Gg\ _ } Ep 
Wair Co _ Re =t- M air oe = Wir em + pe = jap (19) 
Cp Rp Cp Cro Gp Ep 


where C49 is the outdoor air concentration of species A, Gg is the indoor generation rate of 
species A, and E is a new variable — the zone excitation — equal to the sum of species mass inflow 
and generation rates. Defining new constants K2;, i = A, B, C, D equal to the composite of the 
leading constants for each of the well-mixed zone rate expressions above (i.e., as tabulated in 
Table 2) then Equation (19) may be rewritten as: 


C4 — Ky4CaCp C4 E, 
Cp — KopCaCp ajacp\. |) £p 

Wair Cc al M gir ae KC, Cp ne M air Ge 7 Ec (20) 
Cp + KopCaCp Cp Ep 


For a multizone building idealization where all or selected zones are subjected to the same 
bimolecular reaction, the system transport matrix will be the sum of a) the contribution due to zone- 
to-zone flow elements that model dilution dynamics [Kgjj] and b) the contribution due to chemical 
kinetics [Kxjn] which assumes the block-diagonal form, given the chemical kinetics is local to a 
given zone: 

MORAY © ovo 


OW feeIK] Om op 
Kaltes air-2 
[ kin] 0 0 2 : 


0 0 ee M gir-nlK| 


(21) 


for an n zone idealization, if concentration degrees of freedom are ordered by species first and zone 
thereafter as {C47, Cgz, Ccj, Cp1, Ca2, Cp2 ... }?. Here, [K] is the zone kinetics transport 
matrix based, for example, on the second terms of Equation 20 for identical bimolecular kinetics in 
each zone. A similar block-diagonal form will result for termolecular kinetics as well. 


2.1.1.7. The Pseudo Uncoupled Linear Form 
Differential equations governing single-zone chemistry based on the assembly of the kinetic rate 
expressions presented in Table 2 may be written in the following special autonomous form (15) 


that has the form (but not the substance) of an uncoupled, linear system of ordinary differential 
equations — the pseudo uncoupled linear form: 


AC} = {P(Cq, Cpr .--)) — [Car Cos -- KO} 22) 


where {P(C,, Cz, ...)} is a production vector and [L{C,, Cp, ...)| is a diagonal matrix of 
apparent system eigen values (i.e., inverse time constants) identified as the loss rate matrix — both, 
importantly, containing only nonnegative terms. This form is termed autonomous since the 
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production vector and loss rate matrix contain terms assumed to be time independent. Equation 21, 
is a form that, in essence, hides the nonlinearities and the implicit coupling of the system equations 


in the individual terms of the P and L. 


For the system of differential equations used above to model bimolecular kinetics in a single, well- 


mixed zone, Equation 20, the pseudo uncoupled linear form becomes: 


a nal =F yn 0 0 0 
CA “5 a 0 air + KC, 0 0 | cd 
mene M air : Eni io a (23) 
af me KacCaCn + 7g 0 0 = ON ved 
KypC4Cp + iP , uy p Me 


air 


It is important to note that, in spite of the apparent form of this equation: 


e Bimolecular (and termolecular) kinetics couples the dynamic response of the reactants and 


products. 


¢ Bimolecular (and termolecular) kinetics introduce nonlinearities in the resulting system of 
ordinary differential equations, here represented as nonlinear coefficients in a system 
having the form of a uncoupled system of linear equations. 


At any instant, say t, during the dynamic response of the single-zone system being studied, 
however, a given state of reactant and product concentrations will exist 
{C,(t), Cp(t), Ccft), Cp(t)} , thus at that instant in time the behavior of the system may be 
described in terms of a linear, uncoupled system of equations with (instantaneous) species time 
constants (T,4, Tp, Tc, Tp) equal to the inverses of each of the diagonal terms of L: 


1 1 1 
tA = Tw eB Sa ee ee at 
(it + KaaCa(0)} (nit + KoeCa(0)] aa 


providing, of course, the nonlinear production terms for products C and D, KycC,Cp and 


air air 


KypC,4Cpz remain relatively constant (or insignificant) at the instant of time under consideration. 
Thus within the limitations of these caveats and the recognition that all terms are nonnegative, the 
characteristic species time constants of the system are bounded as: 


1 M gir 1 M air Mer 
TA (es > tp = (>; teat = A (24 
: Tae a) ‘B vee a coleaty & a2 eae 


That is to say, the species time constants have an upper bound determined by the dilution time 
constant and, for the reactants, a lower bound that may approach zero when the product of the 
composite rate constants and species concentrations, K24Cg and/or K2pC4, greatly exceed the 
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dilution eigenvalue, Wgj7/Mgir. From a numerical point of view, then, the dilution dynamics will 
tend to place a upper limit on the system time constants (i.e., a lower limit on the systems 


eigenvalue) which will tend to moderate the stiffness of the resulting system equations. 


2.1.1.8. Solution Methods — PSSA and Newton-Raphson Methods 


The nonlinear, coupled, systems of equations that result when bimolecular, termolecular, or other 
nonlinear chemical kinetics are modeled may be solved by a number of techniques. Here we will 
consider two broad classes of methods — methods based on the so-called pseudo-steady-state 
approximation (PSSA) scheme and those based on the Newton-Raphson Method. To make these 
considerations concrete we will apply these methods to the single zone bimolecular example 
considered above and draw general conclusions from these applications. 


PSSA Scheme: The pseudo-steady-state approximation scheme (also identified as the asymptotic 
approximation scheme) is based directly on the pseudo, uncoupled linear form presented above. If 
the production vector P and the loss rate matrix L characterizing the dynamic system’s behavior 
remains constant then an exact solution of the system equations is available (15): 


{C(t +82)} = eHLYCQ)) + [1 - e-SIL][ LP (P) (25) 


This observation has lead to a number of explicit, time-stepping, integration schemes based on 
discrete approximations to Equation 25 (15) such as: 


{Coy} =e PlLAKC,} + [1 — e SIL NLT {P,} (26) 


where a compact notation has been used: {C,} ~ {C(t,)}, [L,] = [L{({C(t,)))] . 

{P,} = {P({C(t,)})} , and t,,,; = t, +t . Given the diagonal form of the loss rate matrix, the 
indicated matrix operations in Equation 26 are computationally trivial, consequently this algorithm 
is computationally very attractive. This scheme is similar to that used by Yamamoto (16) for the 
US EPA/AEERL IAQ model. 


One may, conceivably, cast multizone contaminant dispersal analysis into the pseudo, uncoupled 
linear form (e.g., by including all off-diagonal dilution transport terms in the production vector) 
and achieve computational efficiency. The resulting loss of accuracy and stability when compared 
to more rigorous implicit and explicit approaches may be expected to become a problem however. 


The PSSA scheme has been presumed to provide quick, computationally inexpensive solutions to 
(single-zone) chemical kinetics problems that provide sufficient accuracy for atmospheric urban air 
shed analysis (i.e., where the chemistry is modeled at a very large number of points in the 
atmosphere, over broad geographic regions such as the Los Angeles basin). Verwer and van Loon 
show results, however, that indicate implicit methods utilizing Newton-Raphson integration, while 
more complex, are computationally competitive and appear to be more reliable than the PSSA 
schemes considered (15). 


Newton-Raphson Method: To apply the Newton Raphson Method, we begin by applying a finite 
difference algorithm to the system of differential equations described by Equation 20 to transform 
them to a nonlinear system of algebraic equations that will be solved in a step-by-step manner to 
approximate the dynamic solution of the original system. To this end, Equation 20 may be 


rewritten in autonomous form as: 


Ca \ Ca TRAC ACE E, 
d}) Ca\ _ _ Wair} Ce — KypC Cp 1_} Ep 
dt Vi Ce (“Min Ce : + KycCyaCp ~™M nen wee ie 
Cy Cp + KynCaCp En 
or, in concise notation: 
ce age Warr eG E 
iGo a M.\©} tH {R(C,, Cz)} 8 M4 } (27b) 


where the vectors contained in Equation 27b are defined implicitly by comparison to Equation 27a. 


To transform this system of differential equations, a general semi-implicit scheme will be employed 
based on the following consistent finite difference approximation: 


{Cri} a {C,} wells at) £{C,} a ot £ {Cyst} (28) 


where & is a parameter to control the nature of integration with 0 < o < 1; the time domain has 
been divided into discrete steps t,,, = t, +61; and an abbreviated notation has, again, been used 
as C,, = C(t,,) . Note that this difference approximation corresponds to the Forward Difference 
Scheme for = 0; the Crank-Nicholson scheme for a = 1/2; and the Backward Difference Scheme 
for & = 1. 


Evaluating Equation 27 at time ty, and tn47 and substituting the finite difference approximation, 
Equation 28, we obtain: 


Wai 


Hie }(Cnet} ~ 284 R(Cp4i))) = (1 - (+ 03s) 
(1+ a8:}{ 57 —(E,} + (RUC,})}\ + 08tq7 —{Enat) 


(1 + dt \(Cn} a 


(29) 


This rather complicated expression defines a step-by-step integration algorithm where the vector of 
unknown concentrations at t = tn+7, {C,,,} is approximated from the values of the other listed 
vectors on the right hand side that have known numerical values at time tn4 7 . For our purposes 
here, the sum of these right hand terms will be identified as an effective dynamic excitation 
{Ei} as: 
{E43} = (1 — (1+ 081) 


air 


(Cn) + (14082 Hz (En) + {RUKC,)))} + a8tz1_{E,,1} 30) 


so that Equation 29 will assume the apparently simpler form: 


(31) 


With ao = 0, Equation 31 becomes an explicit integration scheme while with a = 1 it becomes a 
“fully” implicit scheme. For values between these extremes, the algorithm leads to so-called semi- 
implicit schemes. 


The implicit variants of Equation 31 demand the solution of a system of nonlinear equations at each 
time step. This may be achieved using the Newton-Raphson method. To do so, we reformulate 


Equation 31 in residual form as: 


(F((Cns1))} = (1 + 081 


}{Cni2} ~ @BRECp4))} — {Ent} 2) 


Wair 
M air 
and form the Jacobian [f'({C,,,;})] of the residual {f({C,,,,})} with individual terms defined as: 


F fC) = Ares 33) 


and apply the following iterative algorithm based on Taylor’s formula: 
LF Crnea) MCs)? -{Cnea)") + F{Cnes)")} = 0) (34) 


where k is an iteration index. 


To summarize, the proposed algorithm involves a step-by-step integration of the discrete form of 
the system equations (Equation 31) using a Newton-Raphson iterative scheme (Equation 34) at 
each time step to solve the nonlinear algebraic problem defined these semidiscrete equations. To 
implement this algorithm one would follow the steps listed below: 


Given the (approximate) value of the system concentration vector {C,,} at time step n 
fue. t= pa 


1. Form the effective dynamic excitation for the next time step {E41} ; 
2. Initialize the iterative algorithm: 

2.1 Setk=1 

2.2 Set {Cyus}' = {Cy} 
3. Iteratively solve Equation 34 until convergence: 

3.1 Form the current system Jacobian f'({C,,,1)*) . 


3.2 Form the current system residual f({C,,,:}‘) 
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3.3 Solve Equation 34 for ({C,,.;)"*? — {Cu41)*) 


3.4 Update estimate for the system concentration vector: 


: k+1 k k 
{Gun} = ((Cra1) . —{Cnat} )+ {Cast} 
3.5 Evaluate convergence 6C of the system concentration vector by 


evaluating an appropriate norm of the solution from 3.3: 


ae k+1 k 
oC = {C41} } — {Cui} | 
3.6 If SC is not sufficiently small, increment k (i.e., set k = k+/ ) and go to 


3.7 otherwise continue. 


4. Step forward to the next time step: 


4.1 Set {Cri} a Great : 


4.2 Setn = nt. 
4.3 Go to Step 1. 


The system Jacobian [f?] is central to the success of this method. First, the iterative algorithm can 
not proceed unless this matrix is nonsingular. For the system being considered the system 
Jacobian is the sum of a contribution due to the dilution transport and one due to the bimolecular 
kinetics: 

00 —K 44 Cpin41) —K 24 Can) 9 0 

00] _ ose —KpCpin41) —K2BCains1y 9 9 

oe KycCainst) KrcCacneiy 99 
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) W.: 
[R'((Cyst))] = (1 + abr 2) 1 
aur 0 0 
K2pCains1) KapCacnsi 00 


(1 + Ot a) POOL; 4 Cae) Coke 0 0 
81K 25 Can41) (1 + adr ial + 81K > Cains) 0 0 (35) 
: ELOTKGEC pig) = O0iKo- C1 (1 + a6 +f) 0 
—O5tK rp Cgin+1) —O5tK 2p Cains) 0 (1 + adr i) 


At this point little can be said about the general nature of this system Jacobian, although given all 
variables and constants in this Jacobian are nonnegative it may be concluded that this Jacobian will, 
in general, not be diagonally dominant. 


It follows from this example that single-zone systems with any combination of bimolecular and 
termolecular chemistry will result in systems of nonlinear equations that may be solved using the 
Newton-Raphson method and the resulting system Jacobian will be the sum of the dilution 


contribution of the form ( + Cora 


4 [2] and relatively simple additive contributions for each 


air 
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reaction considered (i.e., scaled by adr ). If the individual rate expressions tabulated in Table 2 
are considered to be elemental then the reaction kinetics contributions to the system Jacobian may 
be directly assembled from element Jacobians formed by differentiating the elemental rate 
expressions with respect to each concentration variable involved and evaluating the resulting 
derivatives at the current state of the system concentration. For chemical kinetics described by 
other rate expressions, a solution may be realized using the algorithm outlined above by providing 
specific reaction subroutines that can, for the current state of species concentrations {C,}, “report” 
the rate of production or removal of each species involved (i.e., to form the residual used in 
Equation 34) and the current value of the element Jacobian terms (i.e., to form the system 
Jacobian) to a master solver routine. Furthermore, given homogeneous chemistry is limited to 
each zone these general strategies extends as well to the multizone case. 


2.1.1.9. Initial Concentrations, Outdoor Concentrations, and Steady State Concentrations 


The consideration of indoor homogeneous chemistry will demand, necessarily, the specification of 
initial conditions and outdoor concentrations for all species involved. Outdoor chemistry is 
photolytically driven and as a result is in constant flux, consequently the specification of outdoor 
concentrations will prove difficult for most users and may, in effect, limit the analysis of 
homogeneous chemistry to the research specialists. Furthermore, trace, immeasurable levels of 
some species may have a significant impact on indoor chemistry, yet it may be practically difficult 
to specify outdoor concentrations — again a specialist may be able to devise indirect strategies to 
estimate outdoor levels. If analysis is limited to indoor chemistry that has already been studied 
guidance is available (6, 8, 9, 17, 18) — in fact, code is available (e.g., the RADCAL routines from 
Nazaroff’s program MIAQ4 (19). If, on the other hand, the next generation NIST IAQ Model is to 
be general then general tools should be provided, but what can be done in these cases? 


In atmospheric chemistry many reactions are so rapid — i.e., due to large rate constants and/or large 
concentrations of reactants such as oxygen or nitrogen — that it is often assumed that the reactions 
proceed instantaneously to steady-state. As a result, for such systems of reactions, the 
concentration of some species may be related to the dynamic variation of other species through so- 
called pseudo-steady-state approximations (PSSA) — not to be confused with, but related to, the 
PSSA integration scheme presented above. These approximations assume the form of algebraic 
expressions as, for example, for the photochemical cycle of NO2, NO, and O3 the PSSA 
concentration of the oxygen radical O is: 


or in the atmospheric chemistry of CO and NO, the PSSA concentration of ozone O3 is: 


21 


where, in both cases, M represents a so-called third body such as nitrogen gas. 


These PSSA algebraic expressions may, then, be used to establish the variation of some species in 
terms of other species whose variation may be easier to characterize. That is to say, PSSA 
algebraic approximations may be used to a) relate some initial concentrations to other better known 
concentrations and/or b) relate some outdoor concentrations to other better known concentrations. 
Furthermore, it may prove useful to use PSSA approximations to estimate the time variation of 
some indoor species concentrations to others. Numerically, this would involve mixing algebraic 
with ordinary differential equations demanding other solution strategies than those considered 


above. 


To provide the greatest generality, the next generation CONTAMxx should, therefore, allow the 
user to define algebraic relations between air pollutant species for either a) specification of initial 
concentrations, b) specification of outdoor concentrations, and/or c) possibly, specification of 
steady-state dependencies of some indoor species on others. 


2.1.1.10. Heterogeneous Loss 


Most gases that will be involved in homogeneous chemical reactions within rooms may also be 
expected to interact chemically or physically with building surfaces — that is to say, when modeling 
homogeneous chemistry one should also model the heterogeneous processes that may influence the 
reactants and products considered. For a variety of reasons — that will be discussed in section 2.3 
of this report — it will prove difficult to model the details of these heterogeneous processes 
faithfully, consequently, these processes are commonly modeled as a heterogeneous loss Rad 
(mass-species-A-time~!) using first order approximations known as deposition velocity models: 


Raa = Pai YamdAma|Ca (36) 


where ¥4,,q is the mean deposition velocity (cm-s"! or m-s-!) for species A to surface m , and A,,4 
is the surface area available for deposition of surface m. Note that a sign convention has been 
associated with Rag that is consistent with that introduced for homogeneous gas-phase reaction rate 
Ra such that the mass rate of removal from the zone bulk-air phase is negative. Thus gas-phase 
deposition may be directly added to the zone rate vector {R} discussed above to include deposition 
transport in system equations governing transport of air pollutants in building systems (e.g., 
Equation 20). 


Regrettably, the deposition velocity model is deceptively attractive — one can easily imagine loss to 
surfaces in terms of the “velocity” of a given species (e.g., molecule, radical, or ion) toward the 
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surface. Yet the actual mechanisms involved are likely to be far more complex involving a variety 
of diffusion processes to, through, and within porous surfaces complicated by sorption processes, 
change of phase, and surface chemistry. It should not come as a surprise, therefore, to discover 
that measured deposition velocities show extremely large variation — i.e., the “constant” of the 
deposition velocity model is not observed to be constant at all. Seinfeld (1) reports values of dry 
deposition measured in the outdoor environment that vary over four orders of magnitude for a 
range of chemical species and two orders for a given species yet show a clear trend related to the 
reactivity of the gas species considered. (In atmospheric modeling, dry deposition is distinguished 
from wet deposition, the latter associated with removal of gases due to water droplets in clouds or 
rain. Within the indoor environment wet deposition could conceivably play a role when water 
sprays are present (e.g., in showers) but is otherwise likely not to be important. On the other 
hand, water layers and absorbed water within building materials may be expected to play a key role 


in “dry” deposition.) 


First order deposition velocity models are presently supported by CONTAM93/94 so the question 
of implementation is moot. Providing users with recommended values for deposition velocities is 
another matter. At this point in time it is probably best to simply point them to the literature — 
several investigators have measured and investigated the variation of deposition velocity in the 
indoor environment and have put forward (tentative) theories to explain and predict this variation 
(7, 20-27). 


2.1.1.11. Recommendations for the Next Generation NIST IAQ Model 


At this point there seems to be two directions that could be taken for the next generation NIST IAQ 
Model — either implement specific chemistry —i.e., the indoor homogeneous chemistry that has 
already been studied — or provide tools for the investigation of indoor homogeneous chemistry, in 
general. Nazaroff’s model for indoor NOx gas-phase, homogeneous chemistry is far more 
complete than any other available model and could be used directly to implement the first strategy. 
This is an especially attractive strategy given his model is available as extremely well-documented 
FORTRAN?77 routines within the MIAQ4 program (19). 


Alternatively, general tools could be developed. To comprehensively cover all possibilities, the 
following general types of reactions would have to be considered (with user input specified for 
each type of reaction): 


e Zero Order Reactions: specified in terms of: 


a) a simple rate “constant”, ko, defined to be either a (true) constant or 
temperature-dependent (i.e., as defined by an Arrhenius-type expression or a 
simple linear expression?) given either a schedule of temperature variation or a 


3A simple linear expression of the form k = a + bT may be used to approximate any other temperature dependent 
relationship about some mean operating temperature (i.e., by a Taylor series about that mean). 
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matching temperature (i.e., thermal) solution for the system being considered, 


and 


b) a reactant and product species name(s) (i.e., identification strings), product 
stoichiometry, and reactant and product molecular weights. 


¢ Pseudo-Steady-State Reactions, Equation 9: specified in terms of: 
a) arate “constant” k; defined to be either: 
i) a (true) constant, 


ii) temperature-dependent (i.e., as defined by an Arrhenius-type expression 
or a simple linear expression) given either a schedule of temperature 
variation or a matching temperature (1.e., thermal) solution for the 


system being considered, or 


iii) photochemically dependent, defined in terms of discretely defined 
absorption cross-sections, quantum yields, and actinic irradiance density 
(i.e., Equation 8) input as wavelength dependent schedules by the user. 


b) a “steady-state” concentration [A]s; defined in terms of an actual constant value 
or in terms of a schedule of concentration variation with time (e.g., ambient O2 
concentration in the first case or CQ? in the second case — both available in 
relatively large amounts with the first remaining practically constant and the 
second varying with time of day to local auto traffic conditions, for example), 
and 


c) areactant and product species name(s) (i.e., identification strings), product 
stoichiometry, and reactant and product molecular weights. 


¢ First Order (Unimolecular) Reactions, Equation 1: specified in terms of: 
a) arate “constant” kj defined to be either: 
i) a (true) constant, 


li) temperature-dependent (i.e., as defined by an Arrhenius-type expression 
or a simple linear expression) given either a schedule of temperature 
variation or a matching temperature (i.e., thermal) solution for the 
system being considered, or 


ii) photochemically dependent, defined in terms of discretely defined 
absorption cross-sections, quantum yields, and actinic irradiance density 
(i.e., Equation 8) input as wavelength dependent schedules by the user 
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b) 


a reactant and product species name(s) (i.e., identification strings), product 
stoichiometry, and reactant and product molecular weights. 


Pseudo First Order (Bimolecular) Reactions, Equation 15: specified in terms of: 


a) 


b) 


C) 


a rate “constant” k2 defined to be either a (true) constant or temperature- 
dependent (i.e., as defined by an Arrhenius-type expression or a simple linear 
expression) given either a schedule of temperature variation or a matching 


temperature (i.e., thermal) solution for the system being considered, 


a “steady-state” concentration [B];,; defined in terms of an actual constant value 
or in terms of a schedule of concentration variation with time, 


a selection of either the general reaction type defined by Equation 2 or the 
special case defined by Equation 10, and 


d) reactant and product species names (1.e., identification strings), reaction 


stoichiometry, and reactant and product molecular weights. 


Second Order (Bimolecular) Reactions, Equation 2: specified in terms of: 


a) 


a rate “constant” k2 defined to be either a (true) constant or temperature- 
dependent (i.e., as defined by an Arrhenius-type expression or a simple linear 
expression) given either a schedule of temperature variation or a matching 
temperature (i.e., thermal) solution for the system being considered, 


b) reactant and product species names (i.e., identification strings), reaction 


C) 


stoichiometry, and reactant and product molecular weights, and 


a selection of either the general reaction type defined by Equation 2 or the 
special case defined by Equation 10. 


Third Order (Termolecular) Reactions, Equation 3: specified in terms of: 


a) 


a rate “constant” k3 defined to be either a (true) constant or temperature- 
dependent (i.e., as defined by an Arrhenius-type expression or a simple linear 
expression) given either a schedule of temperature variation or a matching 
temperature (i.e., thermal) solution for the system being considered, 


b) reactant and product species names (i.e., identification strings), reaction 


C) 


stoichiometry, and reactant and product molecular weights, and 


a selection of either the general reaction type defined by Equation 3 or the 
special cases defined by Equations 11 and 12. 


Third Body Combination Reactions, Equation 14: specified in terms of: 
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a) rate “constants” k2, and k2.. defined to be either (true) constants or temperature- 
dependent (i.e., as defined by an Arrhenius-type expression or a simple linear 
expression) given either a schedule of temperature variation or a matching 
temperature (i.e., thermal) solution for the system being considered, 


b) reactant, product, and third body species names (1.e., identification strings), 
reaction stoichiometry, and reactant, product, and third body molecular 


weights, and 


c) specification of a broadening factor F to use. 


¢ User-Specified Reactions: specified in terms of: 


a) a user supplied subroutine that can “report” production or removal rates of 
reactant and species and element Jacobians for the specific reaction for a given 


state of the system concentration vector, and 


b) reactant and product species names (i.e., identification strings), reaction 
stoichiometry, and reactant and product molecular weights. 


c) specification of a broadening factor F to use. 


To provide the greatest flexibility, the user should be allowed to first specify a specific type of 
reaction that may be considered (i.e., in one of the categories presented above), then to indicate in 
which zones each specific reaction is to be applied. This will allow the user, for example, to 
specify different photochemical dissociation reaction for different zones depending on the actinic 
intensity variation within each zone. In addition, as discussed above, it would be most useful to 
provide an algebraic parser that would allow the user to define algebraic relations between gas- 
phase species for either a) specification of initial concentrations, b) specification of outdoor 
concentrations, or c) specification of (some) indoor concentrations in terms of others. 


In addition, to facilitate the consideration of gas-phase homogeneous chemistry the next generation 
of the NIST IAQ Model could provide a) libraries of standard reactions (e.g., NOx chemistry to 
begin with), b) concentrations of components of the standard dry atmosphere, and c) a utility to 
convert between gas-phase water concentrations [H20] and relative humidity RH given schedules 
of or computed values of either (Seinfeld’s Appendix 4.A.3 Calculation of Atmospheric Water 
Vapor Concentration provides one approach (1)). 


2.1.2. Aqueous- & Aerosol-Phase Chemistry 
Aerosol-phase chemistry, especially aqueous aerosol chemistry plays a very important role in the 
atmosphere (see, for example, Chapter 5 of (1)). Recent research points to the possibility that 


aerosol chemistry may play an important role in indoor air quality (7, 11-13, 17, 28, 29). 
Aqueous, acid aerosols may have a more significant impact on health than previously thought and 
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can, if present, have a devastating impact on electronic equipment, building materials, and 
contents. Consequently, we may begin to look to a future need for including aerosol-phase 
chemistry in indoor air quality modeling. 


Two important classes of aerosol-phase chemistry are likely to be important — chemical reactions 
occurring homogeneously within aqueous droplets and chemical reactions occurring on surfaces 
and within the pores of solid particles. Diffusion transport through air-phase boundary layers 
separating aerosol particles from the bulk air-phase play a role in both cases and, for this reason, 
both classes of chemistry may be considered to be heterogeneous (i.e., occurring at the bounding 
surface) relative to the bulk air-phase even though homogeneously distributed through out the bulk 
air phase. Furthermore, gas-liquid absorption equilibrium for the aqueous-phase and gas-solid 
adsorption equilibrium for the solid-phase aerosols will prove important in these classes of 
chemistry — additional characteristics of this class of chemistry that supports the positions that these 
processes ought to be modeled as heterogeneous processes. 


That is to say, while it is premature at this time to propose general modeling tools for indoor 
aerosol-phase chemistry, the modeling tools that will be presented in Section 2.4 for heterogeneous 
processes are likely to be useful for this purpose and should be developed with this application in 
mind. 

2.2. Aerosol Transport & Fractional Particle Filtration 


Professor Bill Nazaroff has presented theory and developed the program MIAQ4 for modeling 
aerosol transport in multizone buildings that has established the state-of-the-art in the field (6, 19, 
22, 23, 30-34). He provides the following general description of the MIAQ4 model: 


“The indoor aerosol dynamics model predicts the time-dependent evolution f particle size 
distributions in buildings. The indoor aerosol is represented by a multicomponent sectional 
description (Gelbard and Seinfeld, 1980). The particle size distribution is divided into a 
number of discrete contiguous sections according to particle diameter. Within each size 
section, the particles comprise one or more classes of chemical constituents. 


The Model employs a flexible multichamber description of a building. ... In simulating the 
evolution of particle size distribution, the model accounts for the influences of direct indoor 
emission, ventilation, filtration, mixing between chambers, coagulation, and deposition 
onto indoor surfaces. Evaporation, condensation, and homogeneous nucleation are not 
included in the model.” (33) 


Nazaroff presents a practically complete description of his aerosol transport model formulation in 

(31). This section will simply follow and review the Nazaroff model formulation presented in this 

paper adding details to aide in the understanding of Nazaroff’s model and modifications needed to 

integrate the Nazaroff model into the framework of the current NIST IAQ model. As such, the 

organization of this section follows that used by Nazaroff —i.e., the level three headings (2.2.1, 
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2.2.2, ...) below are identical to those used by Nazaroff so that both this section and Nazaroff’s 


key paper can be reviewed in parallel. 


Fundamental principles of aerosol science and technology are presented in a number of texts (1, 
35-39). Any of these texts should help to clarify unfamiliar concepts — the texts by Hind and 
Seinfeld are especially useful in this regard. 


2.2.1 Aerosol Representation 


Aerosols are suspensions of solid or liquid particles in air. In general, an aerosol may be expected 
to contain a continuous distribution of particle sizes (i.e., if a statistically significant population of 
particles are present), consequently to characterize an aerosol one must characterize the aerosol’s 
size distribution. To this end, a size distribution function n(Dp) is defined so that the number of 


particles per volume (e.g., number-cm-?) Nj within a particle diameter Dp range “j” or section “j 
AD, = (Dj D,j2) is equal to the integral of the size distribution function over that range: 


Ppj2 
Nj= | nDpaD, 37) 
Pyjl 


The total number of particles per unit volume N is then: 


N = |" nD,)ab, (38) 


In addition to the distribution of number of particles, it is often useful to characterize the particle 
surface area distribution ns, particle volume distribution ny, and particle mass distribution ny. 
Functions for these distributions may be derived from the basic number distribution functions if a) 
it is assumed that all particles are spherical and b) for the mass distribution, all particles have the 
same density Pp (€.g., see Chapter 7 of (1)): 


n(D,) = ™D5n(D,) (39) 
my y= %Dpn(Dp) (40) 
Nm(Dp) = PpgDpn(D,) = Ppny(D,) (41) 


The accumulated volume (fraction) of particles Vj (volume-particles-volume-air~!) in a given 
section j is, then: 


D.; 
v= | nop aD, (42) 
P pil 
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Nazaroff utilizes the mass concentration C4 (mass-particles-volume-air!) of each of several 
(discrete) size fractions as the fundamental independent variable of his model. It then follows from 
the above definitions that the mass concentration of section j or C; is: 


DS; sy ena 
eis Picge oT 48 eee Pie 
C; = if Pps D,n(D,) dD, = IE Pp n,(D,) dD, (43) 
pjl pj 
or, if it assumed that all particles have the same density: 
C; = pV; (44) 


The fundamental independent variable of the NIST IAQ Model is concentration C expressed in 
terms of mass fraction, hence the (mass fraction) concentration Cj (mass-particles-mass-air!) of a 
given section j of particles is related to the (volumetric) mass concentration C; and the volume 


fraction V; of that section as: 


Ct eG MRL y. (45) 


where £,;, is the density of the air. 


Thus the distribution of aerosol particles may be described in terms of the discrete distributions Nj, 
Cj, Cj, and/or V;. While the last three are related by Equation 45 (i.e., assuming constant particle 
density), the link between the discrete number size distributions Nj; and the remaining three 
depends on knowledge of the continuous size distribution function n(Dp) which will, generally, 
not be available. In Nazaroff’s model, however, the distribution of particles within each section is 
assumed to be uniformly distributed with respect to particle diameter or, equivalently, with respect 
to the logarithm of particle mass. If, it is assumed that the distribution function is constant within a 


section rl (D Drp)| = nj ) then: 


D 
pj2 
N; = | n.dD,, = n{Dpj2 - Dp} (46) 


7 
Pi bio a7) 
DiA— D4. 
ao OE TG ( pie ma} zi 2 2 
= Pp 94 Ni Rian eas. 24 N (Dy + Dyno + Don] 
( Pj2 nit} 


and 


4 An outline font will be used to distinguish mass concentration variables C (e.g. mass-particles per unit volume of 
air) from mass fraction variables C (e.g. mass-particles per unit mass of air). 
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In addition to distinguishing concentrations of particles within a given particle diameter section, 
Nazaroff distinguishes chemical components — identified by index k — and chambers or well-mixed 
zones, i, So that Cjz “represents the mass concentration of component k in section j contained 
within chamber i. Furthermore, “components are assumed to be mixed internally, i.e., within a 
section all particles have the same composition.” Finally, Nazaroff places a geometric constraint 
on the width of the section used in computation — the mass of the largest particle within a section 
must be at least twice that of the smallest particle within the section ((31) pages 157-158) or, 
assuming constant particle density: 


3 1/3 ats 
Do 2 2Da Or Dp 2 2°" Dp = 1.2599 Dy (49) 


Although the next generation of CONTAM will base computations on particle concentrations 
expressed in terms of mass fraction Cj (mass-particles-mass-air !), users of CONTAMxx may 
prefer to review results in terms of the discrete distributions of Cj, Vj , or Nj. The relations 
presented in this section define the conversions that will be required to report results in these 
alternative forms. 


2.2.2 [Nazaroff’s] Basic Model Postulate 


Nazaroff presents the form of the system equations that will be assembled from elemental or 
component equations for ventilation, filtration, coagulation, etc. in Equation | of his paper. The 
form used will be recognized as that of the pseudo, uncoupled linear form considered in the 
discussion of homogeneous gas-phase chemistry presented above, Equation 22. This form is the 
basis of the pseudo steady state approximation PSSA (or asymptotic approximation) scheme used 
for approximate, but computationally cheap, numerical integration in analysis of gas-phase 
atmospheric chemistry and is the solution method, in fact, used by Nazaroff (i.e., see Computer 
Implementation Notes on page 161 of Nazaroff’s paper (31)). 


As discussed above, the PSSA scheme is not a rigorous approach to numerical integration and can 
not be expected to provide a reliable, accurate, and/or stable means to solve the (more general) 
types of nonlinear equations that may be expected to be encountered in the next generation of the 
NIST IAQ Model. For this reason, Nazaroff’s solution strategy will not be pursued here. Instead, 
we will extract the elemental and or component mass transport relations presented by Nazaroff and 
reformulate them, as necessary, for assembly in more rigorous solution schemes. 


2.2.3 Ventilation and Filtration 


The advective transport of aerosols in a multizone building may be treated in the same manner as 
that used for gaseous contaminants. That is to say, if it is assumed that within the bulk-air phase 
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of a given zone and at any section along discrete flow paths between zones the concentration of the 
aerosol is uniformly distributed (well-mixed), then the methods to model the advective transport of 
gases presently used may be directly applied to aerosols. Of course, as for advective transport of 
gases, there will be instances where these well-mixed assumptions may be inappropriate — e.g., for 
dense aerosols injected into a zone that will tend to behave more like a cloud that a well-mixed 
suspension of particles or the movement of even trace aerosols through flow passages, such as 
cyclone separators, that tend to separate the aerosol suspension from the gas phase. 


Modeling of the advective transport of aerosols will, however, demand significantly larger systems 
of equations. For example, to faithfully capture the detailed nature of aerosol transport, one may 
have to consider from 10 to 30 sections — Nazaroff used 16 sections in his studies (31) — for each 
component considered. Thus for a single component, the number of equations to be solved may 
be expected to be 10 to 30 times the number of zones! 


From a numerical perspective, it will prove difficult to order these equations to minimize either 
their bandwidth or their skyline and, as a result, sparse matrix equation solving techniques, that 
have fostered the use of implicit solution methods, may not be as effective as they have proven to 
be in simpler contaminant dispersal analysis. Nevertheless, some relatively obvious strategies are 
available to mitigate this problem. From the discussion presented above, it is clear that 
homogeneous chemistry will couple concentration degrees of freedom within zones but not 
between zones. We shall see that sorption transport, source models, and, for aerosols, coagulation 
modeling all have the same characteristic. Furthermore, these transport processes will, most often, 
be characterized by nonlinear couplings. Consequently, it may be most reasonable to order system 
degrees of freedom by zone. If this is done, the contributions of these zone processes [Kzgne ] to 
the system transport matrix /K/J- or, perhaps, more importantly to the system Jacobian — will 
assume a block-diagonal form similar to that described by Equation 21. The patterns of coupling 
between these blocks due to advective (dilution) transport, then, will be identical to the coupling 
that would exist for if a single component. That is to say, the system transport matrix may be seen 
as the sum of dilution transport contributions following a skyline pattern and a block-diagonal form 
associated with the zone processes as: 


[Kail 0 0 [Kail BW [Kepre| 0 0 0 0 
\ [Kail Oe ase 0 ome 0 0 0 
[K] = OO Raul On |Kai|| | 9 9 |Krone| © 0 eh 
0 [Kaill 0 [Kail| -- QO 0 0 0 [Korma ey 
a er oe Kal 0 0 0 ade fret] 


The bandwidth of the zone contribution will approach or equal the number of pollutants considered 
while the advective (dilution) contribution will necessarily have a larger bandwidth (more 
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exaggerated skyline) minimally twice that of the zone contribution but more likely several multiples 
of the zone contribution. Consequently, by ordering the zone numbering strategically, one may 


minimize the overall system bandwidth. 
2.2.3.1 Fractional Aerosol Filtration 


The single component filter model presently available in CONTAM93/94 is a particularly simple 
component. Given a filter of efficiency 1x for air pollutant component k, with a supply and, by 
continuity, an exhaust air flow rate Wgir (mass-air-time-!), the supply and exhaust mass flow rates 
of component k:, wz; and wxq_ respectively (both defined positive for flow into the filter) are 
related to the supply Cy; concentration as: 


Wri = WairCgi (Sla) 


Wo = —W,i fl -—n ) Cki (52a) 


Nk 
K f 
Cis Wair Cie Wair 


Figure 1. Single component filter idealization. 


Expressed in terms of an element equation describing component k mass flow rates between 
discrete spatial nodes associated with the supply and exhaust concentrations, these equations may 


Whi \ _ Wair 0 OF : 
ah r ; Pa —Ny 0 oe (53) 


Demanding the conservation of component k mass, the mass rate of accumulation of component k 


be written as: 


on the filter wzis then: 
Wie = Wail Cri (54a) 


and the total mass of component k accumulated Mj, on the filter during a time interval (ty, to) iS: 


19 : 
Mig = I, Wai Chat (55a) 


These relations are based on an implicit definition of the filter efficiency Nx: 


C.-C 
ie: = aE (56a) 
ki 
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These simple relations may be directly applied to each section j of n sections of an aerosol 


distribution j = 1, 2, ... n by simply adding an addition index for the section designation as: 


Wiki = WairC jxi (51b) 
Wiko = —Wai ll — 1 ig) Cini (52b) 
Wing = War jaC jxi (53b) 
7) 
M ing * Wair jC j iat (S5b) 
1 


and, assuming again that aerosol particles are effectively spherical, of uniform density, and 
uniformly distributed with respect to diameter (i.e., by Equations 45 and 48): 


ee ii Cjto _ Opi Ceo _ Nii Niko _ Vitiz Vito (56b) 
’ C iki C iki N xi V iti 


The next generation of the NIST IAQ model should compute and report the distribution and total 
particulate mass accumulated on filters, Mjx¢ and py M ig given an initial (distribution) of 
accumulated particulate mass input by the user. 


2.2.3.2 Fractional Filtration Efficiency — The Empirical Approach 


In the practical application of the equations presented above a number of issues need to be 
considered — importantly, the variation of filter efficiency with particle size, face velocity, 
component, and /oading,. From results just now being disseminated from a major investigation of 
fractional aerosol filtration sponsored by the US EPA and being carried out by the Research 
Triangle Institute (RTI) (40), data is now becoming available to quantify the first three of these 
issues — the fourth remains problematic. Specifically, RTI is presently releasing data to 
characterize the dependency of the filter efficiency on particle size Dp, air flow expressed in terms 
of face velocity v,;,, and loading — determined by RTI indirectly in terms of the pressure drop 
measured across a filter due to the accumulation of a standard test dust. RTI 1s not attempting to 
discriminate filter efficiency variation with aerosol component type. The RTI data set — a data set 
that will clearly become the definitive (only?) data set available into the near future — can be used 
directly to specify the particle-size and face-velocity dependency of the filter efficiency for one of 
the several filters studied, i.e., n jDpp v,i) for single component analyses and D,j presumed to 
the mean diameter of a given section j. The variation with loading is another matter. 


Ideally, in air quality analysis, it would be best to define loading in terms of the mass of aerosol 
accumulated on the filter. One possibility, would be to define loading as the accumulated sum of 
this mass across all sections 2 M . While this could be computed directly during analysis it is 
not without shortcomings, e. g., the impact of the accumulation of large particles may be expected 
to differ from that of small particles; of solid particles from liquid particles, etc. Nevertheless, 


Bee) 


given the state-of-the-art, this definition of loading is likely to prove practically useful, in that one 
may expect it to be well-correlated with both the change of fractional aerosol efficiency and the 
change in the pressure-drop across the filter. Unfortunately, as presented the RTI data can not be 
directly related to loading defined in this manner — there is some chance, however, that RTI has 
collected but not reported data (i.e., change of weight of filters upon loading) that would make this 


possible. 


Unfortunately, there is an additional complication. The RTI investigators chose to evaluate filter 
efficiency in a somewhat unconventional manner. Instead of applying the conventional definition 
presented above, Equation 56b, they used the following definition: 


C jko 


Jko-background 


lip pa ab (37) 


C ici — © jki-background 

where Cjki-background and Cjko-background ate supply and exhaust “background” concentrations — 
particle concentrations measured when the filters were not subjected to an aerosol-laden supply air 
flow. The measured supply “background” concentrations were always at or near zero, but the 
exhaust “background” concentrations were nonzero for loaded filters. A loaded filter will tend to 
shed particles creating nonzero exhaust concentrations even when supply concentrations are zero, 
thus this definition was employed to “remove the effect of shedding from the computed efficiency” 
((40) page 173). 


Given these less than ideal circumstances, the following recommendations are made for the use of 
measured data for aerosol transport analysis in the next generation of the NIST IAQ Model: 
¢ following Nazaroff’s approach (31), it should be assumed that filter efficiency depends 
not on composition but only particle size — that is, for all components of the same size 
fraction, a given filter should be assumed to have the same filter efficiency, and 
¢ minimally, efficiency data should be defined in terms of user input empirically 
measured profiles of efficiency versus mean particle size n jDzj) , and 
¢ better, users should be able to input profiles for a range of two or more discrete values 
of air face velocity through the filter n jDpp Vain and CONTAMxx would linearly 
interpolate between values, or 
¢ ideally, users should be able to input profiles for a range of discrete values of both air 
face velocity through the filter and loading, i.e., as defined above, n jDy a RD De): ip 
and CONTAMxx would estimate the loading, using Equation 55b, and linearly 
interpolate between values. 


Data may be limited to values at the mean particle diameters of each size range used to represent the 
aerosol — to avoid the computational overhead of interpolation — or, alternatively, at any convenient 


set of mean particle diameters that span the range used to represent the aerosol — thus demanding 
interpolation of values for computation. 


2.2.3.3 Fractional Filtration Efficiency — The Theoretical Approach 


Hinds presents well established theory to predict the fractional filtration efficiency of filters 
composed of fibers and the dependency of this efficiency on face velocity (see (39) Chapter 9). 
The details of this theory will not be presented here — the presentation in Hinds is concise and to 
the point and need not be replicated — but the general features of this theory will be reviewed. 


The fractional efficiency of a fibrous filter may be estimated by considering the theoretical 
efficiency of a single fiber Ey. This single fiber efficiency may be estimated as the sum of five 
contributions: 


contributions due to interception (R), impaction (I), Brownian diffusion (D), interception of 
diffusing particles (DR), and gravitational settling (G). An additional term to account for 
electrostatic attraction may be included, but is often neglected given lack of knowledge regarding 
the state of charge of supply particles. Expression for each of these contributions are available — 
the face-velocity entering the expressions for all but the first contribution — that could be used 
directly to provide an analytical approach to the determination of n jp  Vair)- The results of the 
application of this theoretical approach show the same trends as the measured data reported so far 
by the RTI group for fibrous filters, thus encouraging serious consideration of this approach. 
Hinds also presents related expressions for the pressure drop across fiber filters that define 
pressure-flow relations that could be directly used in the airflow computations of the NIST IAQ 
Model — a very attractive added advantage of this theoretical approach. 


Although the theory is limited to fiber filters — indeed the n {Dz , Voir) Curves reported by RTI for 
other types of filters show quite different trends — and does not account for loading, it would 
appear wise to at least investigate comparisons between RTI measured and predicted profiles of 
n jDyp Vir) using this theory and if this comparison is favorable to implement filter components 
based on this theory for the next generation of the NIST IAQ Model. 


2.2.4 Aerosol Deposition onto [Room] Surfaces 


As noted by Nazaroff, particle deposition to room surfaces is likely to be “of secondary 
importance” from a mass conservation perspective but may be primary when the deposition itself is 
of concern — e.g., for users studying air quality in museums or in facilities with valuable electronic 
equipment. Beyond this concern, however, it is likely that aerosol deposits play a key role in 
sorption transport and the surface chemistry that may be associated with it and their impact on, at 
least, the perceived air quality within buildings (41, 42). 


of) 


The nature of particle deposition is understood, at least, qualitatively to involve the movement of 
aerosol particles through surface or boundary layers separating surfaces from the bulk air phase 
indoors. Rather obviously, the movement of an aerosol particle through surface boundary layers is 
strongly influenced by the detailed nature of room air flows near surfaces. Less obviously, the 
movement of such particles is also influenced by electrostatic charge differences between particles 
and surfaces, thermophoresis, photophoresis, diffusiophoresis, radiation pressure, and Stefan 
Flow (39). The three classes of -phoresis are due to momentum gradients associated with 
temperature gradients in the bulk air phase, caused by radiant absorption of individual aerosol 
particles, and due to concentration gradients, respectively. Stefan flow — the bulk flow of 
evaporating (or condensing) gases away from (or to) surfaces can also transfer momentum to 
aerosol particles accelerating them away (or toward) surfaces. In principle, then, aerosol 
deposition may be expected to be governed by the details of airflow, electrostatic charge, 
temperature, concentration, radiation, and bulk motion of evaporating or condensing gases near 
surfaces — details that must remain largely unknown for practical air quality analysis, given the 
analytical effort that would be required to predict them. Consequently, it has become common 
practice to model the net effect — i.e., the mass rate of deposition of particles in section j of 
component k Rjxq — of these processes simplistically as a first order process (i.e., linearly 
dependent on the bulk air-phase fractional aerosol concentration Cjx) using so-called deposition 
expressions of the form: 


Rita = Paid & ¥méAma)C jk (59) 


where Varad is the mean deposition velocity (cm-s! or m-s-!) for particles in size fraction j to 
surface m , and A,,, is the surface area available for deposition of surface m. Note that a sign 
convention has been associated with Rjxq that is consistent with that introduced for both 
homogeneous gas-phase reaction rate Rg and gas-phase deposition rate R4q , 1.e., for species A, 
such that mass rate of removal from the zone bulk-air phase is negative and rates adding species, 
component, or size section mass to the zone (i.e., production, generation, emission, etc.) are 
positive. 


As noted by Nazaroff, vertical, horizontal upward facing, and horizontal downward surfaces must 
be distinguished. For each of these orientations, the Nazaroff model provides for deposition due 
to: 1) natural convection driven by a boundary layer temperature difference at surfaces, 2) 
homogeneous turbulence in the core of the room, or 3) forced laminar flow over the surface. 
Deposition velocity models for the first and second cases are discussed in the key paper paralleling 
this discussion (31); all three models are derived in earlier papers (22, 30); and comparisons of the 
use of the last two deposition models are presented in (34). 


For forced laminar flow past planar surfaces, Nazaroff used the correlation in his earlier 
publication (30): 


1/2 
Vitd = 0.678D*"( =] a (60) 
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where u,, is the air velocity parallel to the surface beyond the boundary layer, L is the distance 
from the leading edge of the surface, and v is the kinematic viscosity of the air-phase. Later, 
Nazaroff approached the problem by solving boundary layer equations. 


In addition to the relations defining the Nazaroff deposition models, one must integrate these 
relations over available surfaces and within each aerosol size section to compute a surface and 
section average deposition velocity to use in Equation 59. Finally three additional relations are 
needed to use Nazaroff’s deposition models: 


e The coefficient of Brownian diffusivity D can be estimated using the modified Stokes- 
Einstein relation (see (1, 39) page 324): 


idlaCe 


Pp 3TUD, (61) 


where, here, k is the Boltzmann constant (1.381 x 10-23 J.°K-!), T is the absolute 
temperature (°K), is the coefficient of viscosity (e.g., of air), and C, is the 
Cunningham or slip correction factor that is dependent on the mean free path A of the 
air-phase molecules (see (1) page 317): 


C..=-1+4 {1.257 + 0.4exp (- 1?) (62) 


Nazaroff’s function DIFFUS(DP, T) computes D given the particle diameter D, and 
absolute temperature 7. Within this function the mean free path is computed using 
function FREEMP(T) and the so-called Knudsen number Kn = 2/1/ D,, is passed to the 
function SLIP(AKN) which returns the slip correction factor that is based on an 
alternative formulation by Phillips, rather than Equation 62, as noted in the code. 

¢ Nazaroff estimates the turbulence intensity parameter Kz (34) using an empirical 
relation developed by Corner and Pendlebury: 


2 0.073, 0a fmiin > 
K, = 0918 OF om} on 


with ky the Von Kérm4n constant taken by Nazaroff to be 0.4, p is the density of air, u 
is the mean air velocity near a surface, and x, is the direction along a surface in the 
direction of the mean air flow. Here the mean air velocity near a surface is 
indeterminate in meaning and, even if more precisely defined, would be difficult to 
estimate. As a practical compromise for estimating the boundary layer thickness, 
Nazaroff takes x; to be the geometric mean of the dimensions of a surface. 


For the next generation CONTAMxx, two approaches may be taken to implement aerosol 
deposition modeling — the user can provide profiles of (surface and section average) deposition 
velocity (possibly, using plotted results presented by Nazaroff as a guide) and/or the models for 
deposition velocity presented by Nazaroff may be directly implemented. As deposition may be 
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seen as a physical process, it is reasonable to assume that deposition velocity is independent of 
composition. 


Given the current state of the CONTAM94 interface, it would be most reasonable to include data 
associated with deposition with the definition of individual building zones. This would allow 
zone-specific specification of deposition velocity, a level of flexibility that may be expected to be 


required to provide realistic deposition modeling. 


As published, the Nazaroff deposition models are directed to researchers in the field rather than 
professional practitioners, as they contain a large number of relatively esoteric parameters that must 
be specified for use and provide estimates of deposition at specific points along a surface rather 
than a surface-average value. Nazaroff has, nevertheless, implemented these models in MIAQ4 for 
practical application so that it should be possible to use the input requirements and computational 
algorithms implemented within MIAQ4 to add these models to CONTAMxx. One communicates 
with the program MIAQ4 through a command-line interface. The following commands pertain to 
modeling aerosol deposition: 
General 
PRINT DEPOSITION 
This commands the program to report computed average accumulation of aerosols on 
surface due to deposition. 
PRINT FLUX 
This commands the program to report averaged source and sink terms including among 
others deposition fluxes. 
SET WALLS 
This command allows the specification of numbers of surfaces for each zone or 
chamber. 
Dir. Peer: D man, 
SET SURFACE LOSS DEPOSITION 
This command specifies that deposition velocities are to be directly specified. 
DATA DEPOSITION VELOCITY AEROSOL 
This commands allows direct specification of deposition velocities for each aerosol size 
section. 
For. inar Flow Model 
SET SURFACE LOSS FLOW 
This command specifies that the forced laminar flow model is to be used. 
DATA SURFACE FLOW 
This command allows specification of the air velocity outside of the boundary layer 
u., , Surface orientation, and geometry used for the computation of deposition velocities 
based on the forced laminar flow model. 
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Natural Convection Model 

SET SURFACE LOSS CONVECTION 
This command specifies that the natural convection model is to be used. 

DATA SURFACES 
This command allows specification of surface temperatures, surface orientation, and 
geometry used for the computation of deposition velocities for the natural convection 
model. 

DATA HEIGHT 
This command allows specification of chamber (zone) height used in the calculation of 
deposition velocity for the natural convection model. 

DATA THERMOPHORESIS 
This command allows specification of the thermophoresis coefficient for each aerosol 


size section. 
Homogeneous Turbulence Model 


SET SURFACE LOSS TURBULENCE 
This command specifies that the forced laminar flow model is to be used. 

DATA TURBULENCE 
This command allows specification of the turbulence intensity parameter in each zone as 
a function of time. 


The source code of MIAQ4 is extremely well-documented, consequently, it should take little effort 


to “port” this code to the CONTAMxx environment. The following functions implement 
Nazaroff’s model (taken directly from the MIAQ4.DOC file): 


@ OC OO OQ 


DVPE - DEPOSITION: NATURAL CONV., HORIZONTAL, ENCLOSED SURFACE 
DVPI - DEPOSITION: NATURAL CONV., HORIZONTAL, ISOLATED SURFACE 
DVPL - PARTICLE DEPOSITION: FORCED LAMINAR FLOW ALONG SURFACE 
DVPTD - DEPOSITION: HOMOGENEOUS TURBULENCE, DOWNWARD SURFACE 
DVPTU - DEPOSITION: HOMOGENEOUS TURBULENCE, UPWARD SURFACE 
DVPTV - DEPOSITION: HOMOGENEOUS TURBULENCE, VERTICAL SURFACE 
DVPV - DEPOSITION: NATURAL CONVECTION, VERTICAL SURFACE 


Users of these models would specify for each surface: a) the orientation of the surface, b) the 


model to apply to the surface — i.e., forced laminar flow, natural convection, or homogeneous 


turbulence — c) the mean air velocity outside of the boundary layer u.., a value that could be based 
on computed flow results with some imagination, and d) the expected temperature drop across the 


boundary layer of the surface (the zone surface temperature would be specified when defining a 


zone). 


2.2.5 Coagulation 


Aerosol particles in a state of motion tend to coalesce when they collide. This process is known as 


coagulation and may be driven by collisions due to Brownian motion, air velocity gradients (shear) 
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in either laminar or turbulent flow fields, differential gravitational settling, or force fields due to 
van der Waals forces or Coulomb forces. Again Seinfeld presents the fundamental theory related 
to these processes (see Chapter 10 of (1)) and, for practical modeling, we shall use the model 
presented by Nazaroff (6) based on an approach developed by Gelbard and Seinfeld. 


Nazaroff writes: 


“The treatment of aerosol coagulation in the model is best considered in two stages: the 
calculation of the collision frequency between two particles and the integration of these 
probabilities to obtain the growth and loss rates for the component masses within each 


section.” (section Coagulation, page 160 of (6)) 


The details of Nazaroff’s model will not be repeated here — section Coagulation of (6) should be 
reviewed before proceeding. In the Nazaroff model, sectional mean coagulation coefficients "B 
are computed for four basic classes of collisions n = 1, 2, 3, 4 — with two variants, a and D , of 
collision classes 1 and 2 distinguished. These coefficients are then used to compute the 
coagulation (c) rate of mass accumulation Rjjxc of each component k within each size section j of 


zone 1:as: 


p=1q=1 
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where Migir is the mass of air in zone i and all concentrations Cjjz are expressed in mass fraction 
rather than the volumetric mass concentration Cj, employed by Nazaroff (see Equation 45 in 
section 2.2.1 above). 


The structure of this rather complex expression will become more evident if we limit it to the 
analysis of a single component aerosol in a single zone and drop the i and k indices — keep in mind 
that coagulation “transport” is limited to individual building zones anyway: 


40 


EOS) [lc ote |B s|C pC | a 


For this simpler case it is more clearly seen that: 


e the first double summation accounts for the collision of smaller particles — i.e., in size sections 
lower than size section j — that will increase the concentration in the j size section, 


e the next summation within the braces accounts for the collision of particles in the j size section 
with smaller particles that result in loss of particles in section j and a corresponding increase in 
particles in larger sections, 


¢ the third term 4 ? B alc jC; accounts for the loss of particles due to the collision of particles in 
the j size section with other particles in section j and a corresponding increase in particles in 


larger sections, and 


¢ the last term accounts for the collision of particles in the j size section with larger particles that 
result in loss of particles in section j and a corresponding increase in particles in larger sections. 


From a numerical point of view, Equations 64 and 65 both describe processes that depend on 
simple sums of products of zone concentration degrees of freedom — analogous to modeling a 
number of bimolecular homogeneous reactions occurring simultaneously within the zone. Again, 
the rate quantity Rjjx- or Rjc has been defined with the same general sign convention sign 
convention as that used for homogeneous chemistry, gas-phase deposition, and aerosol deposition 
so that the contribution of coagulation to the system transport equations will be similar to that 
presented earlier (e.g., Equation 19). Anticipating considering the use of a Newton-Raphson 
method for the solution of system equations involving coagulation we may define a coagulation 
rate vector for a single zone {R.} = {Ric Roe: R3e --- Ree} ’ for s size sections and a coagulation 
Jacobian [R’-] — with i,j components for the single component case: 

Ri. = i (66) 


J 
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and carefully observing the limits on the summations we obtain: 
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With these coagulation Jacobians in hand, then, system equations governing the dispersal of 
aerosols could be solved using a step-by-step implicit integration scheme with Newton-Raphson 
iteration within each step in a manner similar to that presented above (section 2.1.1.8). More will 
be said about this and other possibilities in Chapter 3 of this report. 


Nazaroff’s model for coagulation appears to be quite general and practically complete (i.e., for 
those cases where evaporation, condensation, and Coulomb attraction may be ignored). 
Furthermore, Nazaroff presents well-documented FORTRAN 77 code in MIAQ4 for the 
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computation of the sectional mean coagulation coefficients "B and for the evaluation of the 
coagulation rate vector {R-} in the following routines: 


C BETCAL - EVALUATE INNER INTEGRAL OF SECTIONAL COAGULATION COEFS. 

C COAGCO - EVALUATE CURRENT COAGULATION COEFFICIENTS BY BILINEAR INT. 
C COAGUL - EVALUATE CONTRIBUTION OF COAGULATION TO LOSS AND GAIN 

Cy OR - COMPUTE SECTIONAL AEROSOL COAGULATION COEFFICIENTS 

C GAUSBT - EVALUATE OUTER INTEGRAL OF SECTIONAL COAGULATION COEFS. 


Therefore, it will be a straightforward task to add this capability to the next generation of 
CONTAM«xx. 


That is not to say, however, that analysis of coagulation transport — or aerosol transport, in general 
— will be computationally inexpensive. It should be clear from the discussion presented above 
relating to aerosol representation, filtration, deposition, and coagulation that the analysis of aerosol 
transport is computationally intensive due in part to the large number of concentration degrees of 
freedom that must be considered and, in part, to the overhead of computing deposition velocities, 
mean sectional coagulation coefficients, and the nested sums associated with coagulation loss, 
gain, and, possibly, Jacobian calculations. If the next generation of the NIST IAQ Model is to 
include realistic aerosol transport analysis, then these computational demands will have to be 
accepted. 

2.3. Heterogeneous Processes 


Heterogeneous processes are, from the perspective of room air quality, mass transport processes 
associated with mass transfer to and from the surfaces bounding the room air. These surfaces may 
be liquid or solid, part of the building construction, furnishings, finishes or decorations, or even 
surfaces of aerosols suspended in the room air. While associated with surfaces, heterogeneous 
processes should not be assumed to be due to simple surface phenomena such as “deposition” or 
condensation to smooth surfaces but rather dependent upon dynamically linked and generally serial 
mass transfer processes including: 


e diffusion transport through boundary layers separating the apparent, exposed surface 
and the bulk air-phase, 


e diffusion transport along the actual surface which, for solids, is likely to be convoluted 
and irregular rather than plane and smooth, 


e diffusion transport within the liquid or solid material below the surface involving, for 
liquids, molecular, Brownian and convective diffusion and, for solids, molecular, 


Knudsen and surface diffusion within the porous interstices, 


¢ sorption transport from the near-surface or pore air-phase — i.e., absorption in both 
solids and liquids as well as adsorption and desorption transport to surface sites within 
the porous structure of solids — and, possibly 
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¢ chemical transformations occurring homogeneously within the solid- or liquid-phase 
below the exposed surface or heterogeneously at surface sites within the porous 


structure of solids. 


Regrettably, the popularity of deposition velocity models have fostered a view of heterogeneous 
transport that is invariably simplistic — the movement of air pollutants from the bulk air-phase to an 
exposed surface implicitly assumed to be smooth, impervious, and of unlimited capacity. These 
models tacitly assume the rate of mass transfer is related to the pollutant concentration in the air- 
phase — an assumption that, we shall see, is not completely unreasonable — and the apparent, 
exposed surface area. For porous solids, the actual surface area available is likely to be several 
orders of magnitude greater than the apparent, exposed surface area and the surface capacity is 
generally limited. Furthermore, the serial nature of mass transfer from the bulk air-phase through 
the surface to the interior below will, generally, depend on several distinct, elemental mass transfer 
mechanisms anyone of which may tend to limit the overall rate of mass transport. 


This section of the report will attempt to look more fundamentally at heterogeneous processes in 
general considering, first, elemental mass transport processes that may play a role and then 
applying models of these elemental processes to develop models to investigate the dynamics of 
sorption transport in rooms and building filtration devices and models to simulate the behavior of a 
variety of pollutant sources. These objectives are somewhat ambitious, given the current state-of- 
the-art, so the present discussion will be limited to simple models of these heterogeneous 
processes. 


2.3.1 Elemental Mass Transport Processes 


To gain a better understanding of heterogeneous transport in general, it is necessary to look 
microscopically at room “surfaces” — at the interface between the room air and the bounding solid 
or liquid, at the air boundary layers separating this interface from the bulk (core) air-phase, and 
deep into the solid or liquid below the interface. To do so, one must necessarily consider these the 
spatial dimension from the bulk air-phase, through the interface, and into the liquid or solid below 
considering the variation of pollutant concentration in both the air-phase and the solid- or liquid- 
phase as one moves along the spatial axis. While in reality these concentrations may vary 
continuously with distance, for our purposes we shall seek to approximate this variation with 
values at discrete locations as indicated in the diagram below. 


Boundary Layer 5 0 x) ao) x3 X4 Xs 


Surface Interface 


Figure 2. Spatial discretization of a region in the vicinity of a surface. 


That is, we shall identify: 


e air-phase concentrations in the room core C,, at the interface of the surface C*, and at 
discrete locations xj i = 1, 2, ... within the solid or liquid below the surface Cj, 
recognizing that for liquids and nonporous solids these air-phase concentrations are 
likely to be negligible while for porous solids they will not, and 


¢ solid-phase (or liquid-phase) concentrations at the interface of the surface C, (or 
C; for liquids), and at the same discrete locations x; i= 1, 2, ... within the solid (or 
liquid) below the surface C,; (or C;; for liquids). 


In the subsequent discussion, we will find it useful to consider two limiting cases, uniform 
distributions of either air-phase, solid-phase, or liquid-phase concentrations — 1.e., 
C,=C,=C3z=..., Cy =Cy=Cyg=... or Cy =Cp=CpR=.... Also, we will most often 
assume the room core concentration is uniform — i.e., the familiar well-mixed zone assumption — 
although this will be done for convenience only, the elemental models presented may be applied as 
well to imperfectly-mixed room models. 


2.3.1.1 Boundary Layer Diffusion e—'—o 


Diffusion through the air-phase boundary layer, sometimes called external diffusion, is intimately 
related to the details of airflow near the surface — details that must remain largely indeterminate due 
to the analytic challenge demanded to predict them. Nevertheless, the contaminant species mass 
transport, w (g-species-s-!), through this boundary layer may be estimated using the discrete 
approximation from boundary layer theory : 


(70) 


where A; is the exposed surface area of the surface (m2), Pair is the density of the air-phase (g- 
air-m-3), and h is the surface-average mass transfer coefficient (m-s-!). For planar solid or liquid 
surfaces within rooms, the exposed surface area would normally be taken as the apparent or 
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projected area of the surface while for airflow over a bed of granules within an activated carbon 
filter, for example, the exposed area would be taken equal to the accumulated surface area of all the 
individual particles — both areas consistent with the correlation used to estimate the surface-average 
mass transfer coefficient. 

Expressed in terms of an element equation describing species mass flow rates between discrete 
spatial nodes associated with the room core and surface air-phase concentrations, these equations 


a = Pair dh ee = \ oa (71) 


where species mass flow from a node into the boundary layer element establishes the positive 


may be written as: 


convention. 


The average mass transfer coefficient may be measured directly, measured indirectly using the 
naphthalene sublimation technique (43), estimated from published heat transfer correlations using 
the so-called heat and mass transfer analogy, or estimated from published mass transfer 
correlations. Correlations are available for airflow through fixed beds of granules (44, 45) that 
may be used in some building sorption filtration systems. 


For the next generation of CONTAMxx, a boundary layer element should be implemented. The 
user would have to specify: 


: the exposed surface area A,, 


° the surface average mass transfer coefficient h or, ideally, would select from a 
choice of correlations for (at least) airflow over a planar surface (43) or through 
fixed beds of granules (44, 45), and 


° the connectivity of the boundary layer element — that is, the discrete nodes 
corresponding to the core and surface concentrations. 


Within the current interface conventions of CONTAM93/94 the specification of the connectivity of 
the boundary layer element and, indeed, the surface node itself presents a small problem. Limiting 
consideration, here, to modeling transport to surfaces of well-mixed zones — i.e., corresponding to 
rooms or, possibly, tanks of a tanks-in-series idealization of an air cleaning devicé (e.g., sorption 
filtration or wet scrubber devices) — then it would be useful to introduce a surface component 
represented by a surface icon. A surface icon could be placed in a CONTAM zone and opened to, 
metaphorically, reveal the micro-structure of the surface — the boundary layer, surface interface, 
and solid or liquid below, similar to Figure 2 above. The core and boundary layer nodes would be 
defined, by default, upon selection and placing of the surface icon in a zone and the user would, 


with the icon open, specify the model parameters for boundary layer diffusion. In addition the 
user would have the option to: 


e define additional spatially discrete nodes corresponding to layers below the surface, 
Cron. Cau, 

¢ define model parameters governing the nature of diffusion into these subsurface layers, 
if any, 

¢ define mass transfer barriers or resistances, 


¢ specify the nature of any equilibrium relations that may exist between air-phase and 
solid- or liquid-phase concentrations at the surface and each of the subsurface nodes, 


and 


¢ specify initial concentrations and/or hypothetical ideal sources at each of the spatially 


discrete nodes. 


That is to say, the introduction of a surface component would allow the user to, in essence, zoom- 
in and define the microscopic mechanisms governing the heterogeneous process of interest in terms 
of the boundary layer process of this subsection — which is fundamental to all heterogeneous 
processes — and the processes presented in the subsequent subsections of this section. This will 
not only allow a great deal of flexibility in modeling individual heterogeneous processes, but will 
allow the user to model more than a single heterogeneous process in each zone. 


2.3.1.2 Mass Transfer Resistance o—n 9 


The use of so-called vapor barriers in building construction is commonplace. These barriers are 
designed to limit the transport of water vapor into building construction and materials and may be 
constructed of thin material layers or applied as paint. They may be expected to also limit the mass 
transport of other species as well and, therefore, may have to be modeled when considering some 
heterogeneous processes. 


Following the example of the boundary layer mass transport, one may model these mass transport 
barriers as simple resistances to mass transfer of negligible thickness and across which the air- 
phase concentration drops from the one concentration C; to another C; as indicated in the diagram 
below: 
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Figure 3. Mass transfer barrier or resistance. 


For an ideal resistance, the contaminant species mass transport, w (g-species-s~!), through this 


layer may be modeled as: 


(72) 


where Ag is the surface area of the barrier (m2), Pgjr is the density of the air-phase (g-air-m-3), and 
Rm is the mass transfer resistance (s-m-!). 


Expressed in terms of an element equation describing species mass flow rates between discrete 
spatial nodes, these equations may be written as: 


Wil ee Os et ees 
fe an pare 1 | C; (73) 


where again species mass flow from a node into the element establishes the positive convention. 


In reality a barrier to mass transfer will have a finite thickness and could, therefore, be modeled as 
the porous solid that it is. This element allows the analyst to avoid considering this level of detail 
and is likely to provide an acceptable level of accuracy for typical vapor barriers. The use of this 
type of element may, however, introduce a small numerical problem. Often, a vapor barrier will be 
located at the exposed surface of a material — at the inner surface of the boundary layer. Asa 
result, one would have to introduce the i and j nodes at this location leaving the room-side node 
without capacitance — that is to say a single algebraic, rather than differential, equation will be 
introduced thereby demanding a solution strategy that can handle mixed algebraic and differential 
equations. 


For the next generation of CONTAMxx, a mass transfer resistance element should be 
implemented. Numerically, it is no different than a boundary layer element so implementation 
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would be trivial. The surface component interface strategy proposed above could be configured to 
allow the user to introduce the required coincident i , j nodes and specify the model parameters. To 
facilitate its use, the user should be allowed to define this element in terms of permeability or other 
parameters commonly used in the vapor barrier field. 


oh wks, Equi LEriRCansirnints O Rieter eatin O 


Sorption, vaporization, and sublimation mass transport are likely to be limited by diffusion 
processes that transport molecular species to or from the actual surface site where the sorption, 
vaporization, or sublimation processes take place. The kinetics of the actual sorption, 
vaporization, or sublimation processes 1s likely to be so rapid that for the purposes of indoor air 
quality analysis they may be treated as instantaneous. These processes are, however, limited by 
equilibrium constraints — that is to say, at the molecular level they may be assumed to be rapidly 
occurring reversible processes characterized by a) an equal exchange rate of molecules across the 
surface interface and, as a result, b) near-surface concentrations that effectively approach the levels 
that would exist at a steady, equilibrium state if the gas-liquid or gas-solid system were closed 
systems. 


Somewhat ironically, therefore, we will not actually describe the mass transport for these key 
heterogeneous processes but instead establish algebraic constraints on near-surface concentrations 
based on empirical equilibrium measurements which, in general, will be identified using the 
functional notation: 


Cy = AC) or Coe AC} (74) 


where f is a function that is unique for each species/solid (or liquid) system. It will prove useful to 
also identify the general inverse relation as well: 


C” = g{Co) or C” = g{Cp) (75) 


As so often is the case, this seemingly benign modeling strategy of representing fast subprocesses 
with pseudo-steady-state, algebraic relations introduces numerical problems. Here the problems 
are threefold. First these types of constraints will result in mixed algebraic and differential system 
equations; second in several important cases — i.e., sorption modeling with effective sorbents — 
these algebraic constraints will introduce numerically challenging nonlinear couplings between 
system degrees of freedom; finally, by themselves, these constraints do not account for mass 
transfer from the air-phase to the solid or liquid-phase and their use may introduce violations of the 
conservation of contaminant species mass. We will attempt to address these problems as we 


proceed. 


49 


Uy soeen OR O 
Gas-phase air contaminants will be soluble in liquids to a greater or lesser degree. If the 
contaminant does not chemically react with the liquid, the absorption of the contaminant may be 


represented as: 


A + liquid  Aeliquid + AH,,) 


where Aeliguid is meant to represent the species A in solution with the liquid-phase and is the heat 
of solution. For dilute solutions, the equilibrium solubility of the gas in the liquid will be governed 
by Henry’s law (see (46) Section 14.12 or (1) Section 5.1): 


[A], = HyDa (76) 


where [A], is the liquid-phase concentration of species A (mass-A-mass-liquid-!), H4 is the 
Henry’s law coefficient with units, here, of (mole-A-1-liquid-!-atm-!), and pa is the vapor pressure 
of species A in the gas-phase above the liquid (atm). Henry’s law coefficients are often related to 
so-called Bunsen absorption coefficients (46). 


Henry’s law coefficients have a temperature dependency defined by van’t Hoff’s relation (see (1) 
for details): 


din(H,) _ AH 
dT  ——- RT? 1) 
which leads to the approximate relation: 
AH alee 
In\H,(T>)) = In\,(7))—-1 = - = 78 
(H4(T2)) = (HAT) (7: | (78) 


If it is assumed that in the gas-phase species A is governed by the ideal gas law, then we can relate 
the gas-phase concentration Cy to the partial pressure as: 
Ma 


(79) 


where Mz is the molecular weight of species A (gram-A-mole-A-!), R is the ideal gas constant 
(8.31441 Pa-m>.°K-!mole-!), and T is the absolute temperature of the gas-phase (CK). The mass- 
fraction liquid-phase concentration of species A is related to the molar concentration as: 


M 
Ci, = Alp (80) 


Thus Henry’s law may be written in a form more directly useful for our purposes as: 


C4 (81) 
or: 
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(83) 


with: 


If the gas contaminant chemically reacts with the liquid in a reversible manner that establishes an 
equilibrium state between the gas-phase species concentration and liquid-phase product 
concentrations it may still be possible to use Henry’s law to describe the equilibrium state using a 
modified Henry’s law coefficient (see Section 5.2 of (1) for examples). 


For the next generation of CONTAM«xx, an absorption equilibrium constraint should be provided. 
The surface component interface strategy proposed above could be configured to allow the user to 
specify the degrees of freedom to be constrained and the model parameters. It would be best to 
also automate the computation of the temperature dependency of this constraint. 


G 
Adso rp tion Oe CQ 


The sorption of a species A to an active site S—* may be represented as (47-50): 


adsorption 


> 
ee 


desorption 


A + S-* AeS—* + AHay 


where A*S—* is the species in a bound state to an active site and AHgq is the heat of adsorption 
released. For closed systems under steady conditions the rate at which adsorbate molecules bind to 
active sites will eventually come to equal the rate at which they are released; the concentration of the 
adsorbate molecules in the air-phase C and the sorbed-phase C; (mass-species/mass-sorbent) will 
remain constant at their respective equilibrium values. These equilibrium values depend on the 
nature of the adsorbate, adsorbent, and the thermodynamic state of the system. For isothermal 
conditions at atmospheric pressure equilibrium relations between these variables — identified as 
adsorption isotherms — may be represented functionally as given in Equation 74. Some isotherm 
models are tabulated below. 


Table 4. Representative adsorption isotherm models. 
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Langmuir C? is the molecular mono- 


layer sorbent concentration. 
K,, is the Langmuir Constant. 
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K per is the BET constant. 


ax* 

Cis the reduced 
concentration: — the ratio of 
the equilibrium to the 
saturation air-phase 


: * 
concentration C,,,,. 


BET 
Freundlich x K,; and n are empirical 
Freundlich coefficients. 
Polanyi DR 9 D is the Dubinin-Radushkevich 
. Gr: parameter. 
Co = C SPA BUN Sere bs 
m iia G C°? is the sorbed concentration 


corresponding to complete 
filling of micropores. 


Modeling adsorption dynamics depends critically on modeling sorption equilibrium accurately. 
This may be achieved only through careful selection of a candidate isotherm model and empirical 
determination of model parameters using carefully measured experimental data. For the sorption of 
trace concentration air-pollutants on building materials the Langmuir model or its low-concentration 
asymptote, the Linear model, are reasonable first candidates. For sorption of water or other 
contaminants with concentrations within one order of magnitude of their saturated values the BET 
model should be considered. While the Freundlich isotherm is often employed for the sorption of 
volatile organic compounds (VOCs) on activated carbon, the Polanyi DR model is the model of 
choice for sorption of VOCs on granulated activated carbon (GAC) where capillary condensation in 
microporous interstices — pore filling —is important. In spite of recent progress in this area (51-54) 
there is very little data available for the sorption of indoor air pollutants, at low concentrations 
commonly found in buildings, on building materials or sorption filtration media. 


The Polanyi DR approach, which is discussed in detail by Shen and Wang (55), is particularly 
attractive because, in principle, it may be used to describe sorption equilibrium for whole classes of 
compounds on a particular adsorbent and the variation of this equilibrium with temperature. This 
relation is most often presented as a so-called characteristic curve of the form: 


W, = Woexp | ~ky(RZ)” n( 2" (84) 


é 


In this form, the (liquid) volume of adsorbate contained in micropores W, (cm3-g-sorbent-!) is 
related to the total volume of micropores available W, (cm3-g-sorbent-!) and the change in free 
energy of the sorbate from the gas-phase to the sorbed-phase as described by the Polanyi potential 
@ = RT \n(C¥c,) , where R is the universal gas constant and T is the absolute temperature of the 
sorbate/sorbent system. Within this Polanyi relation, three constants, Wo, ky (a constant 
associated with the microporosity of the sorbent) and m (which is 2 for activated carbons), relate 
only to the character of the adsorbent and the other terms, the affinity b and the gas-phase 
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concentration at saturation Csq;, relate only to the adsorbate. The affinity provides a relative 
measure of the affinity of the sorbate for the adsorbent and is though to be related to the molecular 
volume Vi (cm3-mol-sorbate"!) of the liquid adsorbate. 


On this basis Shen and Wang have proposed the following general form of the Polanyi isotherm 
relation: 


von mona : 


and have fitted this form to measured data for a set of six hydrocarbons on two sorbents — a 
petroleum-based activated carbon identified as BAC and an organic sorbent based on polystyrene 
(55). These authors provide a comparison between the characteristic curve they developed for the 
BAC activated carbon with several other characteristic curves for other activated carbons. While 
these curves differ, as expected, the curves presented for a variety of coconut shell-based activated 
carbons, commonly used in building applications, are similar to each other and to the curve 
obtained for the BAC. Regrettably, however, all of the curves presented were developed for 
relatively high air-phase concentrations ranges (i.e., on the order of 100 to 20,000 ppm). 


The author adapted the Shen/Wang general formulation by adjusting the sorbate parameters W, and 
K to better fit the few equilibrium measurements available that approach the concentration levels 
found in buildings (i.e., on the order of 1 - 500 ppb) to obtain (56): 


T \2 C.\2 
Cy. = 0.496 pexp { -0.0016(E) in( 2} ] 


é 


(86) 


with Wo = 0.496 cm3-g-sorbent"!, K = 0.0016, and Cs = We p (p is the density of the liquid 
form of the adsorbate, g-cm-3). This general Dubinin-Radushkevich isotherm — plotted for 
isobutyl methyl ketone, toluene, and benzene (with normal boiling points of 117 °C, 111 °C and 
80.3 °C, respectively) — is compared, in Figure 4, with measured data and several (higher level) 
isotherms reported in the literature for room temperature conditions. Figure 5. compares toluene 
isotherms for three temperatures, based on Equation 6, with isotherms reported by Yeh (57). 
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Figure 4. Representative General Dubinin-Radushkevich isotherm compared to: (a) isotherms 
reported in the literature for toluene (solid line) and benzene (dotted line) — (a. on petroleum-based 
activated carbon (58), b. on coal-based activated carbon (59), c. on coconut shell-based activated 

carbon (60), and d. on phenolic-based activated carbon fiber (61)); (b) measured data for 
isobutylmethy] ketone (triangles) and toluene (crosses) on coal-based activated carbon (59) and 
benzene (diamonds) on coconut shell-based activated carbon (62). 
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Figure 5. Representative General Dubinin-Radushkevich isotherm compared to Freundlich 
isotherms reported by Yeh for toluene on petroleum-based activated carbon at 25 °C (a.) and 75 
PC Dye DT 
It is important to stress measured data is simply not available to establish a single sorption isotherm 
for any VOC, or other air pollutant for that matter, at the 0-1000 ppb levels commonly encountered 
in buildings. The strategy used here of adjusting an available, high-concentration level isotherm to 
fit the few low-concentration measured data points available must be considered to be a method of 
last resort and not a recommended procedure for obtaining needed low-concentration range 

isotherms. 


More recently, Wood published correlations for the coefficients of the Polanyi DR isotherm for 
some 140 compounds and 15 different activated carbons (63). While the data is "limited to vapors 
of organic compounds that exists as liquids at ordinary temperatures and pressures" and is based 
on higher concentration level measures than that found in indoor air, this data should prove useful 
in extending the scope of available data to a broader variety of VOCs. 


So-far, the discussion above has been limited to single-component sorption. In practical situations 
building materials and sorption filtration media will be exposed to multi-component mixtures of 
compounds that will tend to compete for available active sites, consequently, multi-component 
sorption equilibrium will have to be accounted for. A variety of multi-component equilibrium 
relations have been developed and their efficacy for industrial applications has been investigated 
(44, 64). Liu has successfully demonstrated the application of a multi-component approach based 
on Polanyi's potential theory and a technique known as component grouping to the sorption of 
low-concentration VOC mixtures on activated carbon (65). 


In general, for multi-component mixtures (under isothermal conditions at atmospheric pressure) the 
equilibrium sorbed concentration for component i , C,, ;, will depend upon not only on the air- 
phase concentration of component i, C, ;, but the air-phase concentrations of all other components, 
as: 

Cr 01 Co 7, Ceo ee} (87) 


The Extended Langmuir Equation developed by Markham and Benton for component molecules 
that do not influence each other provides one example of a multi-component model: 


C,,, ;K7 ;C.; 
Cx 2 as Late (88) 


1+ 2 Ki jC. 


The model parameters are defined for each component as presented above in Table 4. For very 
low air-phase concentrations (i.e., 1 >> 2.K,,C,; ) this multi-component model simplifies to a 
series of n uncoupled single-component Linear isotherms and, therefore, for this limiting case, 
multicomponent competition may be ignored. 


For the next generation of CONTAMxx, an adsorption equilibrium constraint should be provided. 
The surface component interface strategy proposed above could be configured to allow the user to 
specify the degrees of freedom to be constrained and the model parameters. It would be best to 
offer the user a choice of isotherm model to use — e.g., from the table of models presented above — 
and also automate the computation of the temperature dependency of this constraint. 


Vaporization and Sublimation Om Lom) 
Liquid molecules with sufficient kinetic energy can escape the surface of a liquid to enter the gas 
phase and, conversely, gas molecules may enter the liquid phase. For closed systems of a single- 
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component A under steady conditions these reversible opposing processes of vaporization and 


condensation may be represented as: 


A gas ri Alig iD AA yap (89) 


where AH, is the heat of vaporization. Eventually this system will come into a state of 


equilibrium — i.e., saturated conditions in the gas phase will prevail — at which point the saturation 


vap 


vapor may be estimated using the Clausius-Clapeyron equation: 


—AH, 
Pa = | en (90) 


The constant in this expression p,,, is conveniently determined from boiling point data for the 
pure liquid A (see (46) Sections 5.9 and 12.4. (Sublimation of a gas from the single-component 


solid phase follows a similar rule.) 


Assuming, again, ideal gas behavior for A — Equation 79 — we obtain the saturation (equilibrium) 
concentration of species A above the pure liquid A: 


PaoM, _|—-AH 
Ca-sat = 5 RT exp| ere] (91) 


For closed multicomponent systems, Raoult’s law describing the equilibrium behavior of ideal 
solutions may be used to estimate the (partial) vapor pressure p; of each component J: 


Pi = Xipf (92) 


where x; is the mole fraction of component i in the liquid solution: 


x; = CM (93) 
~ C,/M ; 
Assuming ideal gas behavior and combining Equations 90, 92, and 93: 
Pico sehr ivap 
Cg = —— Beg a 
PaiRTL CJM, Re 


where AH;,,, is the heat of vaporization, M; the molecular weight, and Cj is the liquid-phase 
(mass fraction) concentration of component j. For a single-component (pure) liquid Cj; = 1.0 and 
Equation 94 simplifies to 91. 


Note, to apply Equation 94 to multicomponent evaporation modeling, one must track the change of 
liquid concentration Cj; . Guo and Tichenor have put forward a model for emissions from wet 
sources (e.g., VOCs from wet interior architectural coatings) based on the assumption that the 
vapor pressure of a given VOC component will vary in proportion to the mass of the VOC 
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component remaining in the coating at any time (66-68). Their assumption may be considered to 
provide a crude approximation to Equation 94 above. 


D 
2.3.1.4 1D Simple Diffusion @-——® 


Diffusion transport within nonporous solids or liquids may be modeled using Fick's Law with the 
diffusivity D based on molecular or convective transport processes. For one-dimensional diffusion 
into planar nonporous solid sheets or bulk liquids with a single free surface — i.e., as opposed to 
spherical droplets — the resulting partial differential equations assume the form of the classic 1D 
diffusion equation: 


0°C aC 

Nonporous Solids pA,D> i PsAcrs = PAs (95a) 
d°C aC 

Bulk Liquids pA D>> + p/A,r) = pA (95b) 


where Pp, and (, are the density of the solid and liquid respectively (g-m-3), A, and A, are the 
cross-sectional areas of the 1D liquid or solid domain (m2), x is the spatial dimension into the solid 
or liquid (m), and r, and r, are specific generation rates of the contaminant (g-species-g-!) 
introduced to account for the possibility of distributed generation (+) or removal (-) of the diffusing 
species. 


A solution to this partial differential equation may be approximated using Finite Element 
procedures by discretizing the solid or liquid into a number of layers (see Chapter 10 of (69) for an 
introduction to these procedures). Using the simplest 1D finite element approximation based on 
linear shape functions and a lumped capacitance approximation, we obtain element equations that 
describe the species mass transport rate (g-species-s-!) from adjacent nodes, i, j, into the element,, 
w;, and w,, as: : 


dC,; 
wi) _ As8xP 7 Le DAD lial) Cu ~ A,oxP,[ 10 dt (96) 
dt 


where 5x is the element (discrete layer) thickness. To model diffusion transport in liquids we 
simply switch the s subscripts to an / subscripts. 


Higher order finite element approximations and so-called consistent capacitance matrices could be 
used to gain accuracy while maintaining computational efficiency here. It is likely, given the 
current uncertainty in the field regarding species diffusivities and, especially, specific generation 
rates, that this simpler finite element approximation will, however, be sufficient for the time being. 
Numerically, then, subassemblies of elements defined by Equation 96 may be directly assembled 
with equations from existing CONTAM93/94 mass transport elements to form system equations 
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that include this class of diffusion transport. Given these equations are linear first order differential 
equations, CONTAM’s current solution method may be applied to the solution of these equations 
and the only problem that may be expected to arise will relate to the “stiffness” of the resulting 
equations. The ratio of the element capacitance and the element conductance (8x7/ (2D,)} provides 


a nominal measure of the time constant associated with this element. 


To apply these equations to practical problems of analysis, initial conditions must be established 
for all nodal concentrations. Setting these conditions to zero corresponds to modeling an initially 
uncontaminated solid (or liquid). By setting these conditions to an appropriate high (initial) value, 
one may model an emission “source” governed by diffusion transport. For this case, research is 
needed to establish reasonable “appropriate” initial values. For both cases, experimental work is 
needed to establish effective diffusivities and/or to verify available techniques to estimate them. 


2.3.1.5 1D Equilibrium-Constrained Diffusion in Porous Solids e—--—_-e 

Diffusion transport within porous solids invariably involves a complex variety of processes 
including molecular diffusion, Knudsen diffusion, surface diffusion, and, possibly, Poiseuille 
flow (44, 70). Molecular diffusion in the larger pores and Knudsen diffusion in the smaller pores, 
where the mean free path of the diffusing species is limited by the pore dimensions, are diffusion 
processes that occur in series. They may be modeled analogously to Fick's Law through the 
introduction of an effective diffusion coefficient D, (m2-s-!) (see section 4.2.2 of (44) for details). 
Surface diffusion, involving the motion of the diffusing species along the pore wall surfaces, 
occurs in parallel to the molecular and Knudsen diffusion processes. It may, also, be modeled 
analogously to Fick’s law but with the complication that the effective surface diffusivity, in 
general, depends on the local sorbed concentration. Yang asserts: | 


“Surface diffusion can be important, and dominating, in the total flux in porous material, 
provided: (1) surface area is high, and (2) surface concentration is high.” (44) page 113 


For industrial sorption separation processes this is likely to be the case. In indoor air quality 
analysis where air-phase concentrations and material’s specific surface areas are relatively low this 
is not likely to be the case and, therefore, surface diffusion will not be modeled directly. 

Sheet-Like Porous Solids 


Unlike simple diffusion, diffusion in porous solids may be expected to be constrained by the 
equilibrium sorption conditions discussed above in section 2.3.1.3. By assuming a) the sorbed- 
phase concentration remains in equilibrium with the porous air-phase concentration and b) surface 
and solid-phase diffusion are both negligibly small, one may obtain a partial differential equation 
that describes one-dimensional diffusion into sheet-like porous solids, as: 


2 dC, 
pA,D.SS + p,A,r, = pA,ede + PAs (97) 
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where: 
¢ the first term accounts for the diffusion of the air pollutant species in the porous gas- 
phase with p being the gas-phase density (g-air-m-3), 
¢ the second term is included to account for distributed generation of the pollutant species 
as before, 


e the fourth term accounts for the accumulation of the pollutant species in the pore gas- 
phase with € being the porosity of the solid (volume-pores-volume-total-!), and 


e the last term accounts for the accumulation of pollutant species in the sorbed-phase with 
p, being the bulk density of the porous solid. (In other formulations, the solid-phase 
density p; is used so that p, = (1-€)/ep; would replace the bulk density used above.) 


Imposing the equilibrium condition (see Equation 74 above): 


dC 0 
C,= f(C) or S = Bee (98) 


Equation 97 may be transformed to assume the form of the classic 1D diffusion equation: 
a°C Of, \a 
pAD.SS + PeAgrs = [oA.e + pApxe |e (99) 


Again we may apply Finite Element procedures to approximate a solution to this equation by 
discretizing the porous solid sheet into a number of layers. Using the simplest 1D finite element 
approximation based on linear shape functions and a lumped capacitance approximation, we obtain 
element equations that describe the species mass transport rate from adjacent nodes, i, j, into the 
element, w;, and w,, as: 


dC; 
wi\ _ _ AsdxPsrsf1\ , PAsDe[ 1 -1]] Gi ae afl 10 dt 
{wi > ~ 2seietel th 4 Perel i GET, 12 pe + pus} 4 dG tdeann 
cam 


Note that each element has an effective capacity M equal to the sum of the mass of air within the 
discrete layer plus the mass of the porous solid within the discrete layer scaled by the tangent slope 
of the equilibrium relation as: 


A.6x of 
M og = ca +PSDG ti a0 (101) 
For effective adsorbents the gas-phase capacitance will prove negligible: 
of A,oxp, of 
<< pe Mf = oC, 102 
Bel CC ies Oo vim So 7 Clo tS (102) 
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and it becomes numerically strategic to employ an alternative formulation expressed in terms of the 


sorbed-phase concentration, ignoring the gas-phase capacitance altogether: 


dC,j 
We Vera: ALOx Date fu Ma ARUNG WET BCs) + Asef 0]/ “dr \ 93) 
paminane ic fF Sx [-1 1 J) gic,) 2 | On ace 

dt 


The inverse equilibrium model gC si) used in this formulation may be expanded in a first order 
Taylor series g {Cai} =e (aa Fae ,j to “linearize” the resulting equation about the current 


estimates of nodal concentrations and, thereby, effect a solution. 


When the sorption equilibrium condition is defined by the linear model (see Table 4 above), 
Equation 100 (or 103) becomes a system of two linear ordinary differential equations and, again, 
the inclusion of this element equation in the CONTAM framework presents no special problem 
other than numerical “‘stiffness.”’ For all other sorption equilibrium models these equation become 
nonlinear and, consequently, alternate solution methods must be applied. One possibility, possibly 
the most attractive, is again to employ a time-stepping scheme utilizing the Newton-Raphson 
method. If this tack is taken, element Jacobians may be directly evaluated once a particular 
equilibrium sorption model is specified. 


Granulated Porous Solids 


Porous diffusion within the granules or pellets of sorption filtration media may, in a similar 
manner, be modeled using a partial differential equation describing radially symmetric diffusion 
within a spherical sorbent particle. The Linear Driving Force (LDF) model has been developed, 
however, to provide a computationally simpler alternative (44). To form an element equation using 
the LDF lumped-parameter approximation, a single node associated with the surface location of all 
sorption granules and a node associated with the overall mean concentration of all sorption 
granules, C,, must be introduced to yield: 


: dC’, 

w™ S| 1 Fa Ce a dt 
i? sel aM f 104 
= R |-11}\e,f* ™01])\ ac gat 


where Rp is the sorbent particle radius and C,, is the sorbed-phase concentration in equilibrium 


with the near-surface air-phase concentration. An appropriate sorption isotherm may be introduced 
to transform C,, to the near-surface air-phase concentration C’. 


2.3.1.6 1D Convection-Diffusion 

Por I S 

Air pressure differences across sheet-like porous materials will result in airflows through the 
materials and an added mass transport process — convective mass transport. The superposition of 
parallel] 1D convection and diffusion mass transport is governed by the partial differential equation 
given above, Equation 99, with one additional term to account for convective mass transport: 


a2C Of \AC 


PA.De ar + PsAgr, — Wair Ge a loa 1 ayiaape (105) 


where wgir is the mass flow rate of air through the porous sheet (g-air-s"!). Again, Finite Element 
procedures may be used to transform this differential equation into a discrete approximation. 
Using the simplest 1D finite element approximation based on linear shape functions and a lumped 
Capacitance approximation, we obtain element equations that describe the species mass transport 
rate from adjacent nodes, i, j, into the element, w;, and w j 2 aS (see (69) section 10.4): 


dC; 
use) A,OxP51s [1 Cj A ox Of, 10|) dt 
ban Ss {1} + [ear] + [con] cf +2 {Pe + PsacHfo ify ac, f 05) 
dt 
where: 
_ diffusion transport is modeled as: [k aig] = setae (4 nd (106b) 
convection transport is modeled as: Easel = 3E Ki is + de | (106c) 


with the so-called upwind parameter, O<@<1, is introduced to control numerical 
stability/stiffness (at the cost of artificial diffusion) during the solution phase. 


Equation 106 defines a set of linear, ordinary differential equations that may, therefore, be directly 
implemented in the current CONTAM93/94 system. That is not to say, the use of these convection 
diffusion element equations is not without pit falls. In fact, assemblages of these Equations are 
likely to lead to stiff system equations. The upwind parameter may be adjusted, with experience, 
to mitigate these problems, but may introduce error. 


To be complete, Equation 106 would have to be complemented with pressure-flow relations for 
porous sheet-like materials — i.e., to determine wajr. The Poiseuille equation is a good candidate 
for those cases where the flow is likely to be laminar. Manufacturer literature is available for 
porous sheet material used for particle filtration in the turbulent regime. 
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Packed B Porous Adsorbents 


Axial convection-diffusion in packed beds of porous adsorbents (e.g., activated carbon) may be 
modeled with a convection-diffusion equation that includes a term to account for boundary-layer 
mass transport between the bed void-space and the surface of the adsorbent pellets (see (44) 


section 4.1): 


°C aC Abs * aC 
pAeD,, a air. pad a(C-C } alte oer (107) 


where A is the cross sectional area of the bed, € is the bed void fraction, D, is the axial dispersion 
coefficient of the bed, A : is the surface area of the adsorbent pellets per volume of bed, A is the 
average boundary layer mass transfer coefficient for the void-to-pellet interface, and C " is the air- 
phase concentration a the pellet surface. This equation would then be coupled to a model for the 
diffusion transport within the adsorbent pellets — e.g., Equation 104. 

Following the examples presented above, Equation 107 could be modeled with a discrete Finite 
Element convection-diffusion approximation with boundary layer transport elements linking each 


of the nodal, air-phase degrees of freedom to a two-node element corresponding to Equation 104. 
Diagrammatically, the idealization would take the following form: 


Convection 
Diffusion Link 


1 
ar => =i Wi, 
¥ Boundary Layer 


Diffusion Link 
Figure 6. Packed bed adsorber idealization based on an assembly of conduction-diffusion Finite 
Elements, boundary layer elements, and Linear Driving Force elements. 


LDF 
Diffusion Link 


The author has employed a tanks-in-series approach, without recirculation, to model both fixed bed 
adsorbers and distributed sorption filtration devices more commonly used in buildings using 
individual tanks that account for the boundary layer transport to the adsorbent pellets and intrapellet 
diffusion transport modeled with the linear driving force model (56). This alternative tanks-in- 
series approach, diagrammed below, is completely equivalent to the approach suggested by the 
diagram above when full upwinding (?=1) is employed. A tanks-in-series model with 
recirculation between the tanks would provide an equivalent means to model partial upwinding. 
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Boundary Layer 
Diffusion Link 


LDF 
Diffusion Link 


Figure 7. Packed bed adsorber tanks-in-series idealization based on an assembly of filter cells 
models of assemblies of well-mixed chambers with boundary layer elements and Linear Driving 
Force elements. 

To provide the means to model sorption filtration devices, the next generation NIST IAQ Model 
should include the tanks-in-series (family of) model(s). The paper (56) provides all necessary 
details for implementation. This modeling approach provides capabilities equivalent to the mixed 
Finite Element-Lumped Element approach diagrammed above, when configured appropriately, yet 
is intuitively more direct and may accommodate the extended surface, distributed sorption filtration 

devices found more commonly in building applications. 


The tanks-in-series approach leads to linear systems of ordinary differential equations when the 
sorption equilibrium relation upon which the filter cells are based are linear and nonlinear systems 
for nonlinear sorption equilibrium models. Given, effective sorbents are best characterized by 
nonlinear sorption equilibrium relations, realistic sorption filtration modeling will demand 
nonlinear equation solving techniques. 
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2.3.1.7 Surface Chemistry 


Surface chemical transformations are likely to be important for all reactive air pollutants — i.e., 
those air pollutants that may be expected to be involved in homogeneous chemical transformations 
in indoor air. The list of this class of air pollutants is very long and includes nitrogen oxides and 
their radicals, oxygen, ozone, and oxygen radicals, sulfur dioxide and sulfate ions, organic and 
inorganic acids, and a large number of volatile organic compounds. The dynamics of surface 
chemistry is determined not only by the kinetics of the surface chemistry but the diffusion and 
sorption processes that lead reactants to the surface and products away. In the indoor air field 
surface chemistry has most often been modeled using the highly simplified deposition velocity 
approach (see section 2.1.1.10 of this report), even though, more complete models are available 
from the fields of catalysis and surface science. This section will, therefore, present an 
introduction to these models and present element models for the simplest cases. 


Hudson classifies surface reactions into three broad classes (71): 
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© corrosion reactions where the adsorbate reacts with the surface to either form volatile 
products, thereby consuming the surface, or a more or less stable compound — e.g., the 
formation of metal oxides on exposed metals, 

¢ crystal growth reactions where the adsorbate is deposited on the surface extending the 


crystalline structure of the existing surface, and 


¢ catalytic reactions where the surface does not directly react or combine with the 
adsorbate but, instead acts to catalyze a reaction between the adsorbate and either itself 


or another adsorbate. 


In the indoor context corrosion and catalytic growth reactions have, possibly, important 
consequences — the former damaging surfaces, an especially important problem for electronic 
equipment, and the latter altering the chemical composition of the indoor air that may be beneficial 
or detrimental from a health point of view. 


Here we will limit consideration to two relatively simple surface reaction mechanisms — 
unimolecular decomposition and bimolecular interactions — based on the introductions to the area 
given by Castellan and Steinfeld (14, 46). 

Unimolecular Decomposition 

The simplest surface reaction involves the sorption of a single species, say A, to an active site, 
S—* , bonding the species to the site, AeS—* , modifying the bonding potentials of A ‘s electrons 
and thereby altering it’s chemical activity. The bound species may then undergo decomposition to 
form product compounds that may or may not be important from an indoor air quality point of 


view: 


K sl 
AeS—* —> Products 


The rate of decomposition, R4, may be assumed to be proportional to the mass of A sorbed on the 
surface: 


R, = —KyM,C, (108) 


where K,, is the first order rate constant in mass units (i.e., see 2.1.1.5.), Ms is the mass of 
surface material involved, and C; is the sorbed concentration of the adsorbate. The rate of 
generation of product species may be determined from the specific stoichiometry of the surface 
reaction following the same procedures used in section 2.1 .1.5 for homogeneous chemistry. 


If the specie’s sorption is governed by a linear adsorption isotherm C, = K Co , with C” being the 
near-surface air-phase concentration then the rate of decomposition will be first order with respect 
to this near surface air-phase concentration (i.e., assuming isothermal conditions): 


Ry = -(KM,Kp)C” (109) 
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If, in addition, the rate of diffusion transport (e.g., boundary layer and porous diffusion transport) 
to the surface is rapid relative to the rate of reaction then the near surface concentration will 
approach the bulk air-phase concentration, C’ ~ C,, and the rate of decomposition will be first 
order with respect to the bulk air-phase concentration: 


Linear Sorption; Rapid Diffusion R, =- (K MK pjC (110) 


Under these rather restrictive conditions, then, the first order assumption of the deposition velocity 
approach is realized (see Equation 36). When diffusion transport is non-instantaneous, the surface 
chemistry dynamics will be spatially distributed, even though first order relative to near-surface 
concentrations locally and the behavior predicted by the deposition velocity approach will not be 
realized. 


For sorption governed by other equilibrium relations the rate of decomposition will, strictly 
speaking, be nonlinearly dependent on near-surface concentrations, although, for all physically 
consistent equilibrium relations that approach linear behavior at near-zero concentration levels, first 
order behavior will prevail at low concentration levels. For example, the Langmuir isotherm model 
leads to following surface decomposition rate expression: 


* 

Langmuir Sorption Ry = -(KM ee (111) 
L 

which, at low concentration levels K 1Ce << 1 , approaches first order behavior with respect to C” 

and, at high concentration levels K 7Ga >> 1 , approaches zero-order behavior. That is to say, as 

active sites become completely filled — saturated — the rate of reaction approaches a constant rate. 

Even this relatively simple isotherm model leads to rather complex surface decomposition kinetics; 

other models will introduce even more complex behavior. 

Bimolecular Interactions 

A more complex surface reaction involves the sorption of two species, say A and B, to two active 

sites, S—* ,. The bound species may then react to form product compounds that may or may not 

be important from an indoor air quality point of view: 


K 
AeS—* + BeS—x ins Products 


The rate of consumption of, say, species A, R4, may be assumed to be proportional to the product 
of the mass of A, M,C,,, and B , M,C, , sorbed on the surface: 
2 
Rats ias a(MC4\M,Cp,] = -KyM; CasCes (112) 
where K,, is the second order rate constant in mass units. As a bimolecular reaction, this class of 


reaction has all the variant special cases considered in section 2:1 for homogeneous bimolecular 
chemistry. Again, the rate of generation of product species may be determined from the specific 
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stoichiometry of the surface reaction following the same procedures used in section 2.1 .1.5 for 


homogeneous chemistry. 


. . . * : * . 
If both species’ sorption is governed by a linear adsorption isotherm C,=KpC , with C_ being 
the near-surface air-phase concentration then the rate of decomposition will be second order with 
respect to the near surface air-phase concentrations (i.e., assuming isothermal conditions): 


Ryiere (KoM?K 4pK pp\CxCr (113) 


If, in addition, the rate of diffusion transport (e.g., boundary layer and porous diffusion transport) 
to the surface is rapid relative to the rate of reaction then the near surface concentrations will 
approach the bulk air-phase concentrations, C “= C, and the rate of decomposition will be second 


order with respect to the bulk air-phase concentrations: 
Linear Sorption; Rapid Diffusion R, = — ({k 2M 2K ApK iG ACB (114) 


For this bimolecular case, application of the Langmuir isotherm model to each reactant 
independently leads to following surface rate expression: 


C4KarCa \ C8.KarCr 
Independent Langmuir Chemisorption Rx, = -(k 51M sf Piku | cece (111a) 
ALY A BL“B 


which, 


¢ at low concentration levels K XCF <<1&K meges: << 1, approaches second order 
behavior 


¢ at high concentration levels K,4;C,>>1 & Kz;Cpz>>1, approaches zero-order 
behavior, and | 


¢ when one species is at a high concentration level and the other is at a low level 
K4,C4<<1 & Kp,Cp >> 1, assumes a first order behavior. 


Alternatively, accounting for competition of sorption sites yet applying the basic Langmuir 
approach leads to the so-called Langmuir-Hinshelwood model (14): 


oR CuksacGe 
C.K 4rCaK pi Cp : (111b) 


Langmuir -Hinshelwood Chemisorption Ra = —- (k hii 5 
(1 +Ky,C, + Kars | 

The use of these surface chemistry models for practical air quality analysis is certainly marginal at 
the present time since not one relevant surface reaction has been elucidated to the point that the 
mechanism is understood and the rate and sorption equilibrium constants are known — air cleaning 
devices employing TiO? are arguably the exception to this position. On the other hand, it is clear 
that surface chemistry plays a significant role in the transport of chemically active air pollutants 
indoors and is central to the effectiveness of some methods or air cleaning (e.g., removal of active 
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air pollutants using activated carbon or other chemically treated filters). Research results reported 
by Spicer suggest a catalytic conversion of NO2 to NO on some surfaces and results reported by 
Ryan indicate the “deposition” of ozone on latex paints is governed by a bimolecular surface 
reaction with sorbed water (26, 72). These surface chemistry models could then be included in the 
next generation NIST IAQ Model to support research in the area. 


From a numerical point of view, the inclusion of surface chemistry in the NIST IAQ Model will 
lead to systems of nonlinear equations for all but the special unimolecular case defined by 
Equations 108, 109, and 110. Consequently, implementation will require a general purpose 
nonlinear solution procedure. The Newton-Raphson method outlined above presents one 
possibility. 

2.3.2 Applications of the Elemental Transport Equations 


The elemental transport and constraint relations presented above may be assembled — e.g., as 
indicated in Figures 6 and 7 above — in a great variety of ways to model a number of generic 
transport processes. The icons associated with each section title above will be used in this section 
to illustrate a few of the many possibilities. In these icons: 


°@ oO acircle indicates a node - 1.e., a discrete spatial location within a building or 
building component at which pollutant species concentrations will be estimated or, 
more rigorously, with which a local species concentration distribution will be 
approximated (e.g., a node in a well-mixed chamber is associated with the magnitude 
of a uniform distribution of concentration in the chamber; a node of a Finite Element 
approximation based on linear shape functions is associated with the magnitude of a 
concentration distribution that varies linearly from adjacent nodes). 


©@ filled circle indicates a capacitance node — a node that inherently has a capacitance 


associated with it given the defining mass transfer relation. 


-O an open circle indicates zero capacitance node — a node that does not inherently have 


e © a gray circle indicates a dependent, constrained, or specified node — a node 
associated with concentrations that are either dependent or constrained to another node 
or are have values that are specified — from a mathematical point a view a node defining 
concentration boundary conditions rather than system degrees of freedom. 


a solid line indicates a mass transfer link — i.e., mathematical relation that 
establishes the rate of transport of mass from one “node” to another 


© wrens  @ GOtted or gray line indicates a constraint link — 1.e., a mathematical 
relation that establishes a algebraic constraint between linked nodal degrees of freedom 
(concentrations) but, importantly does not directly define a mass transfer relation. 
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Assemblies of the elemental relations presented above define the so-called topology of a given 
mass transport model — a graph of the (sub)system equations associated with the model that 
completely defines (quantitatively) the structure and (qualitatively) the character of the (sub)system 
equations. Given the icon descriptions above, one may anticipate two difficulties that may arise in 
the use of these elemental relations — the practical use of the constraint relations and zero 
capacitance nodes. 

2.3.2.1 Practical Use of Constraint Relations 

The equilibrium constraints presented above do not define mass transport relations by themselves, 
consequently, they can not be directly assembled with the elemental mass transport relations to 
form mass conservation equations governing the system as a whole. Instead, these equilibrium 
constraints allow one to express mass transport relations in terms of alternative concentration 


degrees of freedom. 


For example, consider the simple, single-zone model with a planar sorption sheet shown below: 


Mair R 


air 


Wair C, 


Figure 8. Simple, single-zone model with a planar sorption sheet — diffusion transport assumed to 
be instantaneous (“Equilibrium Adsorption Model’). 

If it is assumed that diffusion transport is practically instantaneous and, therefore, sorbed 

concentrations Cs; remains in equilibrium with the zone air-phase concentrations C — the basis of 

the Equilibrium Adsorption Model (47-49; 73) — then three equations will be associated with this 

model: 


Zone Mass Conservation Wairl + Wags + ee = wi, + R (112) 

Equilibrium Constraint C, = f() or C = g(C,) (113) 
; dC, 

Sorbent Mass Conservation Waas + M5 wi 0 (114) 


where w,q, iS the pollutant adsorption mass transport rate, C, is the outdoor concentration, and R 
is the zone pollutant generation rate. 


From the third equation w,4, = M,@Cs/4; and from the sorption equilibrium constraint relation we 
note: 


dC 
eo Reni aly 
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Consequently: 


_ yp LO ac 


Wads = SI ACs dt (116) 
With this result the first conservation relation, Equation 112, becomes: 
AF(C) 
we Gat ec) + Mgi |e = Waly + R (117) 


Alternatively, using the inverse equilibrium relation, Equation 117 may be rewritten in terms of the 
sorbed-phase concentrations as: 


= W 


see AC steak (118) 


Wair8(Cs) ¥ [, m0 1) ee aC, dt 


That is to say, this equilibrium constraint relations allow the elimination of one conservation 
relation —i.e., it constrains one concentration degree of freedom to another. 


For adsorption and absorption mass transport, we may impose the equilibrium constraint relations 
within the element-assembly framework of the NIST IAQ Model by introducing the following 
mass transport relations: 


Adsorption mass transport from a near-surface air-phase node to an adjacent sorbent node is 
described by Equation 116, reproduced below: 


Wads = yi) O)ac (119) 


Given a solution for the near-surface air-phase concentration C(t) , the sorbed phase concentration 
C(t) may be recovered using the inverse equilibrium relation. 


In a similar manner, absorption mass transport from a near-surface air-phase node to an adjacent 
liquid node is described by: 


Wabs = M,Aee (120) 


where My is the mass of liquid associated with the liquid node and H is the Henry’s law 
coefficient (see Equation 82). Given a solution for the near-surface air-phase concentration C(t) , 
the liquid-phase concentration C,(t) may be recovered using the inverse equilibrium relation. 
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Cc Csat 


Ey r iy Tr r —_—. 
Evaporative mass transport introduces a new complication — loss of mass from a liquid can not be 
treated as a trace mass transport processes. For example, consider the simple, single-zone model 


with an evaporating liquid layer shown below: 


Figure 9. Simple, single-zone model with a evaporating liquid layer — boundary layer diffusion 
transport included. 


If it is assumed that the gas-phase concentration at the liquid surface remains at the saturation value 


C,,, then three equations will be associated with this model: 


Sat 


Zone Mass Conservation Worl Bevan aa Mire = = W,irCo + R (lady 

Boundary Layer Transport Waap = = Paid sA(c -C a (122) 
ane , a(M C}} 

Liquid Mass Conservation > Weyap erat ae (123) 


Where Weyg, is the evaporative mass transport rate, M, is the mass of liquid that varies with time, 


and C; is the liquid concentration of the pollutant species being studied. 
For a pure (single component) liquid the saturation concentration is determine independently of 
concentration degrees of freedom by Equation 91 as: 


PAM Ay, 
Chea x Rt eee wo oz 


and the liquid concentration is unity C; = 1.0. Consequently, we are left with two uncoupled 
equations that govern the behavior of this system: 


WairC + Pairdsh(C— Coa) + Mai SE = WairCo + R- (124) 
and 
dM 
ee = = pir a(C - en (125) 


To generalize: for pure liquid evaporation, the equilibrium condition Equation 91 simply 
establishes a boundary condition for boundary layer mass transport. The system equations would 
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be used to solve for the bulk air-phase concentration and the changing mass of liquid may be 
computed by simple direct (numerical) integration of Equation 125. 


For multicomponent liquid evaporation, we must consider mass conservation of each component i 
as: 


dC, 


air dt ay 


WairC; i Pair sh C;— Car, i + M airCo, i + R; (126) 


where (Strictly) a different boundary layer mass transport coefficients must be evaluated for each 
component and the component saturation concentrations are coupled to all component 


concentrations — i.e., by Equation 94: 


Pico bi ivap 
C.., = eee le 
Hee OY Uae Xu CyM; Rohe! 


Practically, however, one may gain sufficient accuracy by estimating current, ¢,, saturation 
concentrations from liquid component concentrations from a previous integration time step, ft, _ |, 


as: 


Ce sat(lp) = ee ie wee oy oer cage (127) 
PairRTL Cy(ty_/M; RT 
J 


If this is done, then the mass conservation of each liquid component species will remain uncoupled 
as: 

dM); : 

Are Pair Ash(C}—Coaa,) (128) 
where M_; is the mass of liquid component i. Equations 128 may, thus, be directly integrated to 
determine total liquid mass and species concentrations as: 


M, = XM, (129) 
M); 

Graces (alia! 130 

Li M, ( ) 


2.3.2.2 Zero-Capacitance Nodes: 


Certain assemblies may leave some system nodes without capacitance — that is to say nodes whose 
behavior will be defined in terms of algebraic, rather than ordinary differential, equations that will 
be coupled to other system nodes. If such nodes are to be allowed, the next generation version of 
CONTAM will need solution routines for mixed differential and algebraic equations. 


This problem can quite naturally arise at near-surface nodes associated with a boundary layer. For 
example, the filter cell idealization shown above in Figures 7 includes a near-surface air-phase 
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node for which species mass conservation is described in terms of an algebraic equation (see (56) 


for details). 


When mass conservation for an algebraic node is defined by a linear algebraic relation, the node 
may be algebraically eliminated. When a nonlinear algebraic relation is involved this may prove 
difficult. Alternatively, one may formulate models in such a way as to avoid a zero capacitance 


node. 
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3. SOLUTION OF SYSTEM EQUATIONS 


This section will briefly consider methods to solve system equations assembled from the proposed 
new mass transport elements and note special structural characteristics of these system equations 
that may lead to computational strategies that minimize computational effort and storage 
requirements. 

3.1 Solution Methods 


The new mass transport elements introduced in section 2. of this report differ from those currently 
available in the NIST [AQ Model in three critical respects: 


¢ Nonlinearities Several of the proposed elements are based on nonlinear formulations 
(i.e., the mass transport/transformation being modeled depends nonlinearly on system 
concentration degrees of freedom). When assembled to form system equations, the 
resulting system equations will be nonlinear as well. 


e Stiffness: The use of several of the proposed elements will result in system equations 
that may prove to be very stiff. That is to say, at any instant in time the eigen values 
corresponding to the current state of the system equations may have values differing by 
several orders of magnitude. As a result, the garden variety equation solving 
techniques will demand integration step sizes so small as to be impractical — so-called 
stiff equation solvers will be needed. 


¢ Algebraic Equations: As presented, some of the proposed elements may introduce 
algebraic equations that will be coupled to the (more common) ordinary differential 
equations governing the behavior of the systems being studied. 


The earlier generations of CONTAM have used, by numerical standards, first or second generation 
equation solving techniques that are relatively easy to understand and implement. The proposed 
elements will introduce truly significant numerical challenges, consequently the time has arrived to 
turn to the specialist in the field for guidance. A useful and, apparently, up-to-date reference is the 
new edition of J. Stoer and R. Bulirsch’s Introduction to Numerical Analysis: Second Edition 
translated by R. Bartels, W. Gautschi, and C. Witzgall (importantly, of NIST’s Center for Applied 
Mathematics) (74) — a reference that is hardly introductory by most standards. 


Stoer and Bulirsch distinguish the following general types of ordinary differential equations (using 
notation and sign conventions consistent with the systems considered in this report) ordered in 
terms of general difficulty to solve: 


d{C} _ 


Linear Explicit Systems Ratha a [AK{C()} + {FO} (131) 
- d{C} 
Linear Implicit Systems [M ep = -[K|{C()} + {E(} (132) 
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diC 
Linear-implicit Nonlinear Systems [M pecs = {G(t, C(t))} (133) 


Quasilinear-implicit System [M(C yee ce {G(t, C(t))} (134) 


Importantly, the solution of quasilinear-implicit type of system “is still the subject of research,” 
(74) page 495. In addition, mixed systems of differential and algebraic equations may also be 


encountered. 


The current and earlier versions of CONTAM assembles linear element equations to form linear 


implicit system equations of the form shown above: 


HC} 


[MISC = -[KKCW) + {E0) (135) 


where [M] is the contaminant capacitance matrix, [K] is the contaminant concentration matrix, 
{E(t)} is the contaminant excitation vector (e.g., due to internal sources and contaminant 
infiltration), and {C(t)} are the nodal concentration degrees of freedom. For these versions of 
CONTAM, the concentration capacitance matrix has been a diagonal matrix so, in principle, the 
transformation of Equation 135 to the linear explicit would be trivial: 


AC) - [uy [KK co) + (MP {E@) (136) 
The proposed new mass transport elements may be implemented in a number of different ways so 
as to result in system equations of one or another of these general forms. It is, therefore, strategic 
to avoid the more challenging forms. As discussed in section 2.3.2.1 and 2.3.2.2, it should be 
possible to limit element formulations to mass transport relations and to avoid zero capacitance 
nodes so as to avoid algebraic equations. If this is accomplished, yet all other possibilities are 
admitted we would, in general, need to solve system equations of the following form: 


KG = 


[mM (C(I -[K(C@)KKCM} + {E(t, CD)} (137) 


The nonlinear contributions to the capacitance matrix would result from the modeling sorption 
transport to sorbents characterized by nonlinear sorption equilibrium behavior (isotherms). 


By limiting element equations (e.g., for 1D diffusion) to lumped mass formulations and thus 
obtaining diagonal capacitance matrices, we may avoid this difficult class of differential equation 
and address instead the better understood explicit nonlinear form: 


d{C e 
KD = Gaiag(M(CONTIK(CKCO) + [diag(M(COYT HEC, C)} 138) 
For the types of problems that will be considered when using the proposed new mass transport 
elements, however, the system equations may be expected to be very stiff. 
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Stoer and Bulirsch review relatively recent methods to solve stiff differential equations of the 
explicit nonlinear form (technically, their discussion is related to a close cousin of this form the 


autonomous form) and note two important requirements: 


* to assure stability the method must be absolutely stable — as defined mathematically — 


which forces one to choose implicit or semi-implicit methods, 


¢ tO maximize computational efficiency the method must include automatic time stepping. 


Ironically, these authors use one of the chemical kinetics problem addressed above in section 2.1.1 
— the decomposition of ozone to oxygen — as an example of a stiff system. 


Quoting directly from Stoer and Bulirsch, (74) page 491: 


“Taking A-stability into account, one can then develop one-step, multistep, and 
extrapolation methods .... All these methods are implicit or semi-implicit, since only these 
methods have the proper rational stability function. All explicit methods considered earlier 
... cannot be A-stable. The implicit character of all stable methods for solving stiff 
differential equations implies that one has to solve a linear system of equations in each step 
at least once (semi-implicit methods), sometimes even repeatedly, resulting in Newton-type 
iterative methods. In general, the [solution] matrix ... of these linear equations contains the 
Jacobian matrix....” 


The implicit solution method used in the single-zone example of section 2.1.1.8., based on the so- 
called generalized trapezoid rule for integration with Newton-Raphson integration applied at each 
_integration time step provides an example of such a, A-stable implicit method. Importantly, this 
example demonstrates that the system Jacobian may be directly assembled from elemental 
contributions. This method did not include, however, automatic time stepping. 


At this point in time it would seem most appropriate to follow the recommendations of Stoer and 
Bulirsch and investigate the use of one of the methods they proposed. The one-step method 
developed by Kaps and Rentrop “are distinguished by a simple structure, efficiency, and robust 
step control” (74) page 491 appears to be a particularly attractive candidate. As in the example 
presented in section 2.1.1.8, the system Jacobian used in this particular method may be directly 
assembled from elemental contributions. 


Press, Teukolsky, Vetterling, and Flannery, in the second edition of their very useful Numerical 
Recipes in C (75), recommend the Kaps and Rentrop algorithm (identified as one instance of a 
-Rosenbrock Method), using, however, a set of algorithm parameters developed by Shampine. 
Press, et al. provide C routines implementing their preferred variant of the Kaps/Rentrop method 
as well as practical advice in its use (see pages 738-742 of (75)). These authors imply that semi- 
implicit extrapolation methods may be more appropriate for larger systems and/or higher accuracy 
and provide C routines and advice for these methods as well (see pages 742-747 of (75)). 
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Within the routines offered by Press et al., LU decomposition is used to solve the linearized 
system equations within each time step. For large sparse systems, the LU decomposition routines 
(i.e., ludcmp and lubksb) should be replaced by appropriate sparse matrix procedures. These 
authors, again, provide guidance — see section 2.7 of (75). An interesting approach that should be 
considered for larger problems is the biconjugate gradient method presented on pages 84 to 89 of 
(75) — this iterative method is unique in that it may be applied to systems that are not necessarily 


symmetric, positive definite systems>. 


It is difficult to make specific recommendations with great confidence here. Nevertheless, it would 
seem reasonable to consider employing the following C routines provided by Press et al.: 


¢ stiff — the Kaps/Rentrop algorithm using Shampine’s parameters (see pages 738-742 of 
(75)), with the equation solving routines ludcmp and lubksb replaced by: 


e linbcg — the iterative biconjugate gradient method (see pages 84-89 of (75)) to solve the 


linearized system equations within each time step using 


e the indexed storage routines for sparse systems also supplied by Press et al. (see pages 
78-83 of (75)). 


Press et al. have developed their various equation solving routines to be “plug-compatible” to 
facilitate the use of alternative routines. Using this “plug-compatibility” the solution phase may, in 
principle, be directed to one of several solvers depending on the nature of the system assembly. 
For the next generation of CONTAM, the use of nonlinear elements or components may be readily 
monitored so that, in principle, a nonlinear or linear solver may be selected as appropriate. 
Practically, however, the Kaps/Rentrop algorithm may be applied to both linear or nonlinear 
systems — the latter demanding iteration within each time step and the former not. Furthermore, the 
stiffness of a given assembly — linear or nonlinear — will not be readily evident so, again, to be safe 
a stiff solution method such as the Kaps/Rentrop algorithm should be used. Within the 
Kaps/Rentrop algorithm, however, it may be strategic to use a direct (linear) equation solver for 
small systems and the proposed implicit one for larger systems. 


3.2 Structure of System Equations 


Given the diagonal form of the capacitance matrix, the structure of the conductance matrix alone 
becomes relevant when considering the structure of the system equations and its impact on storage 
and solution time requirements. The discussion of section 2.3.3 applies generally to system 
equations assembled from the proposed mass transport elements with a single caveat — (nonlinear) 
equilibrium constrained convection diffusion between building zones will not be admitted. That is 


5 Press et al. note that a variant of this method for symmetric but non-positive definite systems “has been 
generalized in various ways for unsymmetric matrices” but the citations noted have not been investigated. 
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to say, assemblies of any or all the proposed new mass transport elements will couple intra-zone 
degrees of freedom only — providing we do not allow convection diffusion transport through 
adsorbents separating zones so that the system transport matrix can be expressed as the sum of a 
zone-zone dilution transport contribution and an (intra) zone contribution as given in Equation 50 
and reproduced here: 


[Kai] 9 0 [Kail]: 0 beg 0 0 0 0 
0 [Kg] 0 0 ~~ 0 Oeri| KF] oeotaw ily 0 0 
[K] Ji [tO} GMO [Kaill Outs [Kail] i 0 0 [esa] 0 0 
0 [Kail 9 [Kail © Gal AON Goda kag <1 0 
eal MEO a 06 Or 0 cal Ke 


providing, system degrees of freedom are numbered in the order discussed in section 2.3.3. 


The dilution contribution will be linear while the zone contribution will include nonlinearities if any 
of the nonlinear element types are employed. That is to say, system nonlinearities will be isolated 
to the zone matrices of the block-diagonal zone conductance contribution and will, therefore, have 
a relatively limited bandwidth. Furthermore, updating the system Jacobian will involve only these 
block-diagonal zone contributions and will, in general, be limited to just a few degrees of freedom 
within each individual zone conductance matrix. 
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4. USER INTERFACE STRATEGIES 


To implement the proposed new mass transport elements and components it will be necessary to 
devise additional user interface conventions to complement those presently established by the 
current generation of the CONTAM program —- CONTAM93/94. Currently, to communicate input 
CONTAM93/94 uses a graphic user interface called a sketchpad to, in essence, define the 
topology ®of a given building idealization and data entry screens to associate data with individual 
icons of the sketchpad diagram. In many instances a cascade of several data screens are associated 
with a given icon of the sketchpad diagram that are organized in a tree-like hierarchy with the root 
data screen directly associated with the icon. In some instances, the data entry screens support 
graphical input of data or graphical representation of input data (e.g., for operation schedules). To 
communicate computed results to the user, CONTAM again provides both a graphical 
representation of results on a sketchpad diagram of the building and cascades of data screens that 
provide both tabular data and plotted results. 


Additional interface conventions will need to be developed for each of the three general areas 
considered in this report — homogeneous chemistry, aerosol transport and fractional particle 
filtration, and heterogeneous processes. 


4.1 Homogeneous Chemistry 


Section 2.1.1.11 of this report outlines the basic input data that would be required to define 
homogeneous chemical reactions within building zones. Two reasonable alternatives should be 
considered for the specification of this chemistry — definition at the project level or at the zone level 
of data input. 


4.1.1 Project-Level Definition - 


The general reaction rate expressions summarized in Table 2 above are expressed in terms of the 
product of a) the mass of air in a given zone Mj, b) a rate “constant” Ko, K), or K2, and c) the 
product of zero of more (current) zone concentrations. Furthermore, for realistic modeling a given 
homogeneous chemical reaction should be expected to occur in each and every zone of a building 
idealization. Consequently, reactants, products, stoichiometry, and rate “constants” could be 
defined at the project level for all zones for each of the reactions to be modeled. Rates of 
transformation would then be computed specifically for each zone given Mir. 


The specification of rate “constants” introduces a few complications. Temperature-dependent rate 
“constants” could also be defined generally at the project level and evaluated for each specific zone 
during computation given either specified zone temperature schedules or computed zone 


© Topology here is the geometric configuration of a given assembly of elements and/or components used to idealize 
a given building/HVAC system including the relative position of these elements or components and their 
connectivity (e.g., as defined by discrete airflow paths). 
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temperatures (1.e., for a future generation version of CONTAM that models the coupled thermal- 
airflow problem). Photolytic rate constants depend on discretely defined distributions of 
absorption cross-sections, quantum yields, and actinic irradiance density. The first two 
distributions may, again, be input at the project level for all zones in general. The actinic irradiance 
density depends on the illumination level and spectral distribution — which may be expected to vary 
from zone to zone and with time — thus, this should be specified uniquely for each zone in terms of 


time schedules. 


Specified and dependent concentrations add another wrinkle. As discussed in Section 2.1.1.9, it 


would be most convenient to allow: 


a) the specification of specific outdoor concentration time histories (i.e., boundary 
conditions for the independent concentration degrees of freedom) as constants or 
discrete time histories (e.g., for O2, N2, CO2, 03, NO, NOdg, etc.), 

b) the specification of other outdoor concentrations in terms of these boundary conditions 
(e.g., the specification of algebraic relations resulting from PSSA approximations), 


b) the specification of specific indoor concentration time histories as constants or discrete 
time histories (e.g., for H2O), and/or 

c) the specification of dependent indoor concentration in terms of computed independent 
zone concentrations (e.g., the specification of algebraic relations resulting from PSSA 
approximations). 


To achieve these objectives, a distinction should be made between independent and dependent 
concentration degrees of freedom at the project level. That is, when specifying the contaminants to 
be included in a given analysis, the user should be asked to distinguish independent and dependent 
concentration degrees of freedom. The specification of outdoor (independent) concentrations, 
presently provided for in CONTAM93/94, should be extended to allow the specification of 
dependent outdoor concentrations in terms of these independent concentrations. Finally, the user 
should be able to specify either dependent concentration time histories or equations defining the 
dependency of dependent concentrations on the computed (independent) zone concentrations. 


To maintain a graphical link between icons on the sketchpad and the associated data screens icons 
should be provided for the following: 

An icon associated with the general, project-level specification of reactions to be 
considered. This icon has a double box to distinguish it from icons associated with a 
specific location in the sketchpad diagram. Reasonably, project-level icons could be 
arranged along one edge (e.g., the top edge) of the sketchpad window to graphically 
distinguish them from the building diagram. 

o An icon associated with the specification of time schedules of actinic irradiance 
density discrete distributions. 
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[c] An icon associated with the specification of independent concentration time histories 
or their dependency on dependent concentration degrees of freedom. This icon would 
be automatically placed in the outdoor “zone” waiting for user definition and could be 


placed by the user in (selected) building zones as required. 


A representative sketchpad diagram including these icons is shown below: 


The icon at the upper left indicates homogeneous chemistry has been defined. The O icons in 
selected zones indicate actinic irradiance density has been defined for these zones; its absence in 
some zones indicates the zero actinic irradiance. The El icons in indoor and the outdoor zones 


provides the link to data screens used to specify dependent concentrations. 


In most instances dependent concentration relations (e.g., PSSA approximations) may be expected 
to defined identically in all zones. For this reason, it may be strategic to allow definition of these 
relations globally — i.e., at the project-level, perhaps through the introduction of a [te] icon — 
leaving the [Cl icon to be associated only with specification of concentration time histories for both 
dependent and, possibly, independent concentrations in selected zones. 


“Selecting” the icon would “open” a root data screen for reaction definitions. This root screen 
would reasonably provide a summary of the reactions that have been defined allowing on the order 
of 25 to 50 distinct elementary reactions and would allow either selection from a list of standard 
elementary reactions (TJ) or definition of new reactions (-!) as suggested by the data screen shown 
below: 
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Homogeneous Chemistry 


React. 1:[<define new>] fev@e React. 21° roomed 
React. 2:[ <define new>] TL .1 React. 22: hee 
React. 3:[<define new>] TL 1 React. 23: alge 


React, 4;[<define new>] ML J React, 24: Mie 


etc. etc. 


If the user selects to define a new reaction he/she would then enter a cascade of data screens 


organized as outlined in Section 2.1.1.11 of this report. 


4.1.2 Zone-Level Definition 


As presented, homogeneous chemistry is limited to well-mixed zones, consequently this chemistry 
could be defined uniquely for each zone. A single icon placed within a zone could indicate the 
modeling of this chemistry within a given zone. While this approach allows greater flexibility, it 
will generally require unnecessary duplication of input and contribute to the congestion of icons on 
sketchpad diagrams. Alternatively, homogeneous reactions could be defined globally as in Section 
4.4.1 and turned on or off selectively through the use of the zone icon. 


4.2 Aerosol Transport and Fractional Particle Filtration 


The interface for aerosol transport will demand first the specification of the aerosol representation — 
the number of size fraction bins and their particle size limits (i.e., within the constraint of Equation 
49). Reasonably this specification should be done globally and could be associated with a project- 
level icon, say, a . Here the users should be able to select to represent the distribution in terms 
of Nj, Cj, Cj, or Vj (see Section 2.2.1). 


4.2.1 Fractional Particle Filtration 


As presented in Section 2.2.3, two approaches to modeling fractional particle efficiencies are 
available — an empirical approach based on measured fractional particle efficiencies of filters and a 
theoretical model that is appropriate for fibrous filters. Not discussed in this report, are theoretical 
models for particle deposition in narrow flow paths that effectively act as filters. As presently 
implemented, CONTAM93/94 allows the user to specify single filter efficiency values to all 
discrete flow paths. For those cases where an aerosol is specified globally, then reasonably 
fractional particle efficiencies would have to be specified for all filters. The addition of duct 
modeling tools in the next generation of CONTAM provides another means to specify particle 
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filtration. To avoid redundancy and to improve the correspondence between the actual 
building/HVAC system and the CONTAM idealization, I would propose the following: 
° A filter icon, available in the set of icons used to define duct systems, that, when 
selected, allows direct specification of the distribution of filter efficiencies or the 
specification of theoretical model parameters for the fibrous filter model. 


° A flow path icon shaded to indicate it acts as a filter that, when selected, allows direct 
specification of the distribution of filter efficiencies or the specification of theoretical 
model parameters for the particle deposition model. 


Direct specification of fractional particle efficiencies should be supported by either (or both) 
graphical representation of tabular input or direct graphical input (1.e., by “clicking” data on a 
graph). Following the standard convention, this data should be represented on a semi-log plot of 
efficiency versus particle diameter in the range of 0.1 um to 10 Um as shown below: 
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For the near future a single distribution as shown may be sufficient. For more refined analysis, 
distributions at different filter face velocities and at different loadings may be necessary with 
CONTAM linearly interpolating between values automatically during computation. 


4.2.2 Aerosol Deposition 


Aerosol deposition is most reasonably defined uniquely for each zone. This could be achieve by 
either using the current zone icon or a new deposition icon, say, CF to access appropriate data 
screens. The former choice avoids screen clutter, the latter makes the modeling of deposition 
evident. In either case, the associated data screens should allow either the direct specification of 
deposition velocities for each size bin or the implementation of Nazaroff’s deposition models. 
Surface areas would have to be associated with each deposition distribution, if the former path is 
taken, or, for the Nazaroff models, the orientation of the surface models would have to be 
specified. 
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4.2.3 Coagulation 


Coagulation is an homogeneous process, not unlike homogeneous chemistry, that is described 
generally for all zones and made specific to a zone by the scaling parameter of the mass of the air in 
the zone Mir (i.e., see Equation 64). Consequently, the user need do no more than specify the 


aerosol representation to effect coagulation modeling. 


4.2.4 Aerosol Source Models 


Aerosols may be generated within a zone due to in-zone occupant activity, due to resuspension 
resulting from cleaning, for example, or due to tracking of dust from one zone to another. Naive 
empirical models have been proposed for these sources that have not been considered in this report 
— possibly these models could be used here. Alternatively, the user may simply specify time 
histories of aerosol generation rates, by fraction sizes, in much the same way that (molecular) 


sources are now specified. 


4.3 Heterogeneous Processes 


The heterogeneous processes discussed in Section 2.3 of this report are unique in that spatial 
discretization at the microscopic level is associated with the models of the processes considered. 
Furthermore, heterogeneous transport is inherently zone-specific. For these reasons, it may be 
reasonable to provide an icon, say | that provides the user with the means to access a micro- 
sketchpad — as if zooming in on a zone surface to view the microscopic character of the surface. 
At the micro-sketchpad level the user may be offered a choice of standard heterogeneous models — 
€.g., evaporation or sublimation from a single component reservoir, one of three Axley sorption 
transport models, etc. — or the ability to define a unique 1D spatial discretization similar to those 
represented in Figures 2, 3, 6, 7, 8, and 9 of this report. General finite element analysis programs 
provide models for accomplishing the latter. 


For the modeling of gas-phase air cleaning devices based on microscopic models of heterogeneous 
transport processes, it would be best for most users to be able to simply select from a library of 
device components. CONTAM, would therefore, assume the responsibility to automatically 
choose an appropriate spatial discretization and provide the user with appropriate response results 
(i.e., supply and exhaust concentrations of the device rather than internal, spatially distributed 
concentration results). 


As the approach to modeling heterogeneous processes is unusual, if not innovative, the application 
of the proposed elemental models to various heterogeneous processes will have to be more 
exploratory and, consequently, the development of interfaces to support these process models will 


have to follow these explorations. 
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